Files

479 lines
22 KiB
Matlab

function DataTable = nappBackground(DataTable,varargin)
%NAPPBACKGROUND Removes background from photoelectron or XAS spectra.
% This function determines and removes the background from PE or XAS
% spectra using mean subtraction or curve fitting (linear or exponential).
% Background correction can be performed on kinetic energy, binding energy,
% or X-ray photon energy scales. Results are optionally normalized by
% dividing instead of subtracting.
%
% SYNTAX
% ------
% DataTable = NAPPBACKGROUND(DataTable, energyStart, energyEnd)
%
% DataTable = NAPPBACKGROUND(DataTable, energyStart, energyEnd, 'Name', Value)
%
% EXAMPLES
% --------
% % Simple mean subtraction with kinetic energy scale
% DataTable = nappBackground(DataTable, 368, 370);
%
% % Linear fit on selected rows with binding energy scale
% DataTable = nappBackground(DataTable, 368, 370, ...
% 'tableSelection', tmpDataSelection, ...
% 'energyScale', 'binding', ...
% 'mode', 'linfit', ...
% 'showFigures', 'off');
%
% INPUT ARGUMENTS
% ---------------
% DataTable table
% Data table containing spectra and energy scales. Must contain
% one of the following variables: 'spectrum_counts',
% 'spectrum_countsPerSecond', or 'spectrum_countsNorm'.
%
% energyStart — Start energies for excluded ranges
% Numeric vector. Each element pairs with energyEnd at the same position.
% The positions have fixed roles:
%
% Position 1 (required):
% Lower bound of the main excluded domain (e.g. the peak).
%
% Position 2 (optional): ALWAYS the edge cutoffs, never a peak region.
% energyStart(2): trim the low-value spectrum edge — excludes all
% data points below this value.
% Use NaN if no low-end trim is needed.
% If you want positions 3+ but no edge cutoff, pass NaN here.
%
% Positions 3+ (optional):
% Each pair [energyStart(k), energyEnd(k)] defines an additional
% arbitrary excluded domain (e.g. a secondary peak or artefact).
%
% energyEnd — End energies for excluded ranges
% Numeric vector, same length as energyStart (or padded with NaN).
%
% Position 1 (required):
% Upper bound of the main excluded domain.
%
% Position 2 (optional): ALWAYS the edge cutoffs.
% energyEnd(2): trim the high-value spectrum edge — excludes all
% data points above this value.
% Use NaN if no high-end trim is needed.
%
% Positions 3+: upper bounds of additional excluded domains.
%
% Note: energyStart(2) and energyEnd(2) refer to values on the sorted
% x-axis (prepareCurveData always returns xData ascending), so the
% convention is identical for kinetic and binding energy scales.
%
% Examples:
% % One peak only
% nappBackground(..., 109, 116, ...)
%
% % One peak + low-end cutoff only
% nappBackground(..., [109 106], [116 NaN], ...)
% nappBackground(..., [109 106], [116], ...) % legacy shorthand, equivalent
%
% % One peak + high-end cutoff only
% nappBackground(..., [109 NaN], [116 119], ...)
% nappBackground(..., [109], [116 119], ...) % legacy shorthand, equivalent
%
% % One peak + both edge cutoffs
% nappBackground(..., [109 106], [116 119], ...)
%
% % Two peaks + both edge cutoffs
% nappBackground(..., [109 106 107], [116 119 108], ...)
%
% % Two peaks + low-end cutoff only (NaN = no high-end trim)
% nappBackground(..., [109 106 107], [116 NaN 108], ...)
%
% % Two peaks, no edge cutoffs (NaN NaN at position 2)
% nappBackground(..., [109 NaN 107], [116 NaN 108], ...)
%
%
% NAME-VALUE PAIR ARGUMENTS (optional)
% ------------------------------------
% 'tableSelection' logical array (default: all rows)
% Logical index to filter rows in DataTable.
% Example: tableSelection = DataTable.quality > 0.8
%
% 'tableVariable' char vector (default: 'spectrum_countsPerSecond')
% Spectral variable to process. Must be one of:
% - 'spectrum_counts'
% - 'spectrum_countsPerSecond'
% - 'spectrum_countsNorm'
%
% 'mode' char vector (default: 'mean')
% Background subtraction method:
% - 'mean': subtract mean intensity in energy range
% - 'linfit': subtract linear regression (poly1 fit)
% - 'exponential': subtract exponential fit (exp2)
% - 'mean, divide': divide by mean (normalization)
% - 'linfit, divide': divide by linear fit
% - 'exponential, divide': divide by exponential fit
%
% 'energyScale' char vector (default: 'kinetic')
% Energy axis to use for background calculation:
% - 'kinetic': uses scientaEkin_eV column
% - 'binding': uses bindingEnergy_eV column
% - 'xRay': uses xRayEkin_eV column (for NEXAFS data)
%
% 'showFigures' char vector (default: 'on')
% Display fit visualization:
% - 'on': show fit plot with data and excluded points
% - 'off': suppress figure output
%
% OUTPUT ARGUMENTS
% ----------------
% DataTable table
% Updated data table with background-corrected spectra. The table
% variable specified by 'tableVariable' is overwritten with
% corrected data. If mode contains 'divide', results are written
% to 'spectrum_countsNorm' instead.
%
% Three log columns are added/updated:
% - logBackgroundMode: fitting method used
% - logBackgroundRange: [energyStart, energyEnd] for each row
% - logBackgroundDate: timestamp of processing
%
% NOTES
% -----
% - Empty or missing spectra in tableVariable will trigger an error.
%
% SEE ALSO
% --------
% inputParser, fit, fittype, prepareCurveData, excludedata
% Define conditions and defaults of input
defaultTableVariable = 'spectrum_countsPerSecond';
expectedTableVariable = {'spectrum_counts','spectrum_countsPerSecond','spectrum_countsNorm'};
defaultMode = 'mean';
expectedMode = {'mean','linfit','exponential','mean, divide','linfit, divide','exponential, divide','interpolate'};
defaultTableSelection = true(size(DataTable,1),1);
defaultShowFigures = 'on';
expectedShowFigures = {'on','off'};
defaultEnergyScale = 'kinetic';
expectedEnergyScale = {'kinetic','binding','xRay'};
% Parse input
p = inputParser;
addRequired(p,'DataTable',@istable);
addRequired(p,'energyStart',@isnumeric)
addRequired(p,'energyEnd',@isnumeric)
addOptional(p,'tableSelection',defaultTableSelection,@islogical);
addParameter(p,'tableVariable',defaultTableVariable,...
@(x) any(validatestring(x,expectedTableVariable)));
addParameter(p,'mode',defaultMode,...
@(x) any(validatestring(x,expectedMode)));
addParameter(p,'showFigures',defaultShowFigures,...
@(x) any(validatestring(x,expectedShowFigures)));
addParameter(p,'energyScale',defaultEnergyScale,...
@(x) any(validatestring(x,expectedEnergyScale)));
parse(p,DataTable,varargin{:});
% Pad shorter exclusion vector with NaN so both always have equal length.
% This ensures backward compatibility when only one vector has 2 elements.
energyStart = [p.Results.energyStart, NaN(1, max(0, numel(p.Results.energyEnd) - numel(p.Results.energyStart)))];
energyEnd = [p.Results.energyEnd, NaN(1, max(0, numel(p.Results.energyStart) - numel(p.Results.energyEnd)))];
% check input (all elements of input variable need to have non-zero size)
if ~all(cell2mat(cellfun(@size,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable),'UniformOutput',false)),'all')
disp('The table variable appears to be empty. I quit now')
error('The table variable appears to be empty.')
end
% add empty table variables (if they don't exist)
existColumn = strcmp('lognappBackgroundRange',DataTable.Properties.VariableNames);
if ~any(existColumn)
DataTable.lognappBackgroundRange = NaN(size(DataTable,1),2);
end
existColumn = strcmp('lognappBackgroundMode',DataTable.Properties.VariableNames);
if ~any(existColumn)
DataTable.lognappBackgroundMode = cellstr(repmat(' ',size(p.Results.DataTable,1),1));
end
% add empty variable to table
existColumn = strcmp('spectrum_countsNorm',DataTable.Properties.VariableNames);
if ~any(existColumn)
DataTable.spectrum_countsNorm = cell(size(DataTable,1),1);
end
% repeat this for all input.
for ii = 1:size(p.Results.DataTable(p.Results.tableSelection,:),1)
if startsWith(p.Results.mode,'linfit')
ft = fittype( 'poly1' );
elseif contains(p.Results.mode,'exponential')
ft = fittype( 'exp2' );
elseif contains(p.Results.mode,'interpolate')
ft = 'linearinterp';
opts = fitoptions( 'Method', 'LinearInterpolant' );
opts.ExtrapolationMethod = 'linear';
opts.Normalize = 'on';
end
if contains(p.Results.energyScale,'kinetic')
if ~startsWith(p.Results.mode,'mean')
[xData, yData] = prepareCurveData(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii},p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
% Convention for excludedPoints: true = point is KEPT for fitting.
% excludedata(...,'domain',[a,b]) returns true for points OUTSIDE [a,b].
%
% Fixed layout of energyStart / energyEnd vectors:
% Position 1 : main exclusion region (e.g. the peak) — required
% Position 2 : edge cutoffs — NaN on either side = no cutoff there
% energyStart(2): exclude xData < energyStart(2)
% energyEnd(2) : exclude xData > energyEnd(2)
% Use NaN if you have positions 3+ but no edge cutoff.
% Positions 3+ : additional arbitrary exclusion regions
% Position 1: main exclusion domain
excludedPoints = excludedata(xData, yData, 'domain', [energyStart(1), energyEnd(1)]);
% Position 2: edge cutoffs (NaN = skip that side)
if numel(energyStart) >= 2 && ~isnan(energyStart(2))
excludedPoints = excludedPoints & ...
excludedata(xData, yData, 'domain', [xData(1), energyStart(2)]);
end
if numel(energyEnd) >= 2 && ~isnan(energyEnd(2))
excludedPoints = excludedPoints & ...
excludedata(xData, yData, 'domain', [energyEnd(2), xData(end)]);
end
% Positions 3+: additional arbitrary exclusion regions
nPairs = min(numel(energyStart), numel(energyEnd));
for k = 3:nPairs
if ~isnan(energyStart(k)) && ~isnan(energyEnd(k))
excludedPoints = excludedPoints & ...
excludedata(xData, yData, 'domain', [energyStart(k), energyEnd(k)]);
end
end
if startsWith(p.Results.mode,'interpolate')
opts.Exclude = ~excludedPoints;
else
opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', ~excludedPoints );
opts.Display = 'Off';
end
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
if contains(p.Results.showFigures,'on')
% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
plot( fitresult, xData, yData, ~excludedPoints );
% legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'kinetic energy [eV]', 'Interpreter', 'none' );
ylabel( 'PE intensity', 'Interpreter', 'none' );
grid on
title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
end
if endsWith(p.Results.mode,'divide')
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}./fitresult(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
else
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii});
end
elseif startsWith(p.Results.mode,'mean')
b_1 = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= energyStart,1);
b_2 = find(p.Results.DataTable(p.Results.tableSelection,:).scientaEkin_eV{ii} >= energyEnd,1);
if b_2 == 1
b_2=10;
end
b_mean = mean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(b_1:b_2),'omitnan');
tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false); %#ok<*AGROW>
end
elseif contains(p.Results.energyScale,'binding')
if ~startsWith(p.Results.mode,'mean')
[xData, yData] = prepareCurveData(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii},p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
% Set up fittype and options.
% Convention for excludedPoints: true = point is EXCLUDED from fitting.
% ~excludedata(...,'domain',[a,b]) returns true for points INSIDE [a,b].
%
% Fixed layout of energyStart / energyEnd vectors (same as kinetic):
% Position 1 : main exclusion region (e.g. the peak) — required
% Position 2 : edge cutoffs — NaN on either side = no cutoff there
% energyStart(2): exclude xData < energyStart(2)
% energyEnd(2) : exclude xData > energyEnd(2)
% Use NaN if you have positions 3+ but no edge cutoff.
% Positions 3+ : additional arbitrary exclusion regions
% prepareCurveData sorts xData ascending, so xData(1) is always the
% numerically smallest value regardless of kinetic/binding axis direction.
% Position 1: main exclusion domain
excludedPoints = ~excludedata(xData, yData, 'domain', [energyStart(1), energyEnd(1)]);
% Position 2: edge cutoffs (NaN = skip that side)
if numel(energyStart) >= 2 && ~isnan(energyStart(2))
excludedPoints = excludedPoints | ...
~excludedata(xData, yData, 'domain', [xData(1), energyStart(2)]);
end
if numel(energyEnd) >= 2 && ~isnan(energyEnd(2))
excludedPoints = excludedPoints | ...
~excludedata(xData, yData, 'domain', [energyEnd(2), xData(end)]);
end
% Positions 3+: additional arbitrary exclusion regions
nPairs = min(numel(energyStart), numel(energyEnd));
for k = 3:nPairs
if ~isnan(energyStart(k)) && ~isnan(energyEnd(k))
excludedPoints = excludedPoints | ...
~excludedata(xData, yData, 'domain', [energyStart(k), energyEnd(k)]);
end
end
if startsWith(p.Results.mode,'interpolate')
opts.Exclude = excludedPoints;
else
opts = fitoptions( 'Method', 'NonlinearLeastSquares','Exclude', excludedPoints );
opts.Display = 'Off';
end
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, opts );
if contains(p.Results.showFigures,'on')
% Plot fit with data.
figure( 'Name', 'untitled fit 1' );
plot( fitresult, xData, yData, excludedPoints );
% legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit 1', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'binding energy [eV]', 'Interpreter', 'none' );
ylabel( 'PE intensity', 'Interpreter', 'none' );
grid on
title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
ax = gca;
ax.XDir = 'reverse';
end
if endsWith(p.Results.mode,'divide')
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}./fitresult(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii});
else
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii});
end
elseif startsWith(p.Results.mode,'mean')
b_1 = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= energyStart,1);
b_2 = find(p.Results.DataTable(p.Results.tableSelection,:).bindingEnergy_eV{ii} <= energyEnd,1);
if b_2 == 1
b_2=10;
end
b_mean = mean(p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}(b_2:b_1),'omitnan');
tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false);
end
elseif contains(p.Results.energyScale,'xRay')
if contains(p.Results.mode,'linfit') % LINEAR FIT
% Fit: 'untitled fit 1'.
[xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
% Set up fittype and options.
ft = fittype( 'poly1' );
excludedPoints = excludedata( xData, yData, 'domain',[energyStart energyEnd] );
fo = fitoptions('Exclude', excludedPoints);
% Fit model to data.
[fitresult, ~] = fit( xData, yData, ft, fo );
if contains(p.Results.showFigures,'on')
% Plot fit with data.
figure( 'Name', 'fit pre-edge' );
plot( fitresult, xData, yData, excludedPoints );
xlabel('Photon Energy eV')
ylabel('partial A-M.E.Y.')
grid on
title(p.Results.DataTable(p.Results.tableSelection,:).name{ii},'Interpreter','none')
end
tmp{ii,1} = p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii}-fitresult(p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii});
elseif contains(p.Results.mode,'mean') % mean
if p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(1) < p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(end)
% prepare Data
[xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
% MEAN
b_1 = find(xData >= energyStart,1);
b_2 = find(xData >= energyEnd,1);
b_mean = mean(yData(b_1:b_2));
elseif p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(1) > p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}(end)
% prepare Data
[xData, yData] = prepareCurveData( p.Results.DataTable(p.Results.tableSelection,:).xRayEkin_eV{ii}, ...
p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable){ii});
% MEAN
b_1 = find(xData <= energyStart,1);
b_2 = find(xData <= energyEnd,1);
b_mean = mean(yData(b_2:b_1));
else
disp('Warning: Issue with NEXAFS data or energies')
end
tmp(ii,1) = cellfun(@(x) x-b_mean,p.Results.DataTable(p.Results.tableSelection,:).(p.Results.tableVariable)(ii),'UniformOutput',false);
end
end
end
% fill selected rows
if endsWith(p.Results.mode,'divide')
DataTable(p.Results.tableSelection,:).spectrum_countsNorm = tmp;
else
DataTable(p.Results.tableSelection,:).(p.Results.tableVariable) = tmp;
end
% write log to table
DataTable.lognappBackgroundMode(p.Results.tableSelection,:) = cellstr(repmat(p.Results.mode,size(p.Results.DataTable(p.Results.tableSelection,:),1),1));
DataTable.lognappBackgroundRange(p.Results.tableSelection,:) = repmat([energyStart(1) energyEnd(1)],size(p.Results.DataTable(p.Results.tableSelection,:),1),1);
%--- Define variables for the Universal Logging Block ---
numRowsTable = height(DataTable);
processIdx = find(p.Results.tableSelection);
% -------------------------------------------------------------------
% 1. Dynamically get this function's name
thisFuncName = mfilename;
colDate = sprintf('log%sDate', thisFuncName);
colCaller = sprintf('log%sCaller', thisFuncName);
% 2. Probe for column existence, create empty string arrays if missing
if ~ismember(colDate, DataTable.Properties.VariableNames)
DataTable.(colDate) = repmat("", numRowsTable, 1);
end
if ~ismember(colCaller, DataTable.Properties.VariableNames)
DataTable.(colCaller) = repmat("", numRowsTable, 1);
end
% 3. Get current time as a string
currentTimeStr = string(datetime('now', 'Format', 'yyyy-MM-dd HH:mm:ss'));
% 4. Determine the calling script using the call stack
stackInfo = dbstack;
if length(stackInfo) > 1
callerScript = string(stackInfo(2).file); % The script that called this function
else
callerScript = " "; % Called directly from the command window
end
% 5. Write the log data only to the rows being processed
DataTable.(colDate)(processIdx) = currentTimeStr;
DataTable.(colCaller)(processIdx) = callerScript;
% -------------------------------------------------------------------------
end