Code covered by the BSD License

# Background correction

### Vincent Mazet (view profile)

30 Apr 2010 (Updated )

Compute an estimation of the background (aka baseline) in chemical spectra.

backcor(n,y,ord,s,fct)
```function [z,a,it,ord,s,fct] = backcor(n,y,ord,s,fct)

% BACKCOR   Background estimation by minimizing a non-quadratic cost function.
%
%   [EST,COEFS,IT] = BACKCOR(N,Y,ORDER,THRESHOLD,FUNCTION) computes and estimation EST
%   of the background (aka. baseline) in a spectroscopic signal Y with wavelength N.
%   The background is estimated by a polynomial with order ORDER using a cost-function
%   FUNCTION with parameter THRESHOLD. FUNCTION can have the four following values:
%       'sh'  - symmetric Huber function :  f(x) = { x^2  if abs(x) < THRESHOLD,
%                                                  { 2*THRESHOLD*abs(x)-THRESHOLD^2  otherwise.
%       'ah'  - asymmetric Huber function :  f(x) = { x^2  if x < THRESHOLD,
%                                                   { 2*THRESHOLD*x-THRESHOLD^2  otherwise.
%       'stq' - symmetric truncated quadratic :  f(x) = { x^2  if abs(x) < THRESHOLD,
%                                                       { THRESHOLD^2  otherwise.
%       'atq' - asymmetric truncated quadratic :  f(x) = { x^2  if x < THRESHOLD,
%                                                        { THRESHOLD^2  otherwise.
%   COEFS returns the ORDER+1 vector of the estimated polynomial coefficients
%   (computed with n sorted and bounded in [-1,1] and y bounded in [0,1]).
%   IT returns the number of iterations.
%
%   [EST,COEFS,IT] = BACKCOR(N,Y) does the same, but run a graphical user interface
%   to help setting ORDER, THRESHOLD and FCT.
%
% - V. Mazet, C. Carteret, D. Brie, J. Idier, B. Humbert. Chemom. Intell. Lab. Syst. 76 (2), 2005.
% - V. Mazet, D. Brie, J. Idier. Proceedings of EUSIPCO, pp. 305-308, 2004.
% - V. Mazet. PhD Thesis, University Henri PoincarĂ© Nancy 1, 2005.
%
% 22-June-2004, Revised 19-June-2006, Revised 30-April-2010,
% Revised 12-November-2012 (thanks E.H.M. Ferreira!)
% Comments and questions to: vincent.mazet@unistra.fr.

% Check arguments
if nargin < 2, error('backcor:NotEnoughInputArguments','Not enough input arguments'); end;
if nargin < 5, [z,a,it,ord,s,fct] = backcorgui(n,y); return; end; % delete this line if you do not need GUI
if ~isequal(fct,'sh') && ~isequal(fct,'ah') && ~isequal(fct,'stq') && ~isequal(fct,'atq'),
error('backcor:UnknownFunction','Unknown function.');
end;

% Rescaling
N = length(n);
[n,i] = sort(n);
y = y(i);
maxy = max(y);
dely = (maxy-min(y))/2;
n = 2 * (n(:)-n(N)) / (n(N)-n(1)) + 1;
y = (y(:)-maxy)/dely + 1;

% Vandermonde matrix
p = 0:ord;
T = repmat(n,1,ord+1) .^ repmat(p,N,1);
Tinv = pinv(T'*T) * T';

% Initialisation (least-squares estimation)
a = Tinv*y;
z = T*a;

% Other variables
alpha = 0.99 * 1/2;     % Scale parameter alpha
it = 0;                 % Iteration number
zp = ones(N,1);         % Previous estimation

% LEGEND
while sum((z-zp).^2)/sum(zp.^2) > 1e-9,

it = it + 1;        % Iteration number
zp = z;             % Previous estimation
res = y - z;        % Residual

% Estimate d
if isequal(fct,'sh'),
d = (res*(2*alpha-1)) .* (abs(res)<s) + (-alpha*2*s-res) .* (res<=-s) + (alpha*2*s-res) .* (res>=s);
elseif isequal(fct,'ah'),
d = (res*(2*alpha-1)) .* (res<s) + (alpha*2*s-res) .* (res>=s);
elseif isequal(fct,'stq'),
d = (res*(2*alpha-1)) .* (abs(res)<s) - res .* (abs(res)>=s);
elseif isequal(fct,'atq'),
d = (res*(2*alpha-1)) .* (res<s) - res .* (res>=s);
end;

% Estimate z
a = Tinv * (y+d);   % Polynomial coefficients a
z = T*a;            % Polynomial

end;

% Rescaling
[~,j] = sort(i);
z = (z(j)-1)*dely + maxy;

a(1) = a(1)-1;
a = a*dely;% + maxy;

end

% delete lines below if you do not need GUI

function [z,a,it,ord,s,fct] = backcorgui(n,y)

% BACKCORGUI   Graphical User Interface for background estimation.

% Initialization
z = [];
a = [];
it = [];
ord = [];
s = [];
fct = [];

order = 4;
threshold = 0.01;
costfunction = 'atq';

% Main window
hwin = figure('Visible','off','Position',[0 0 750 400],'NumberTitle','off','Name','Background Correction',...
bgclr = get(hwin,'Color');

% Axes
haxes = axes('Units','pixels');

% Buttons OK & Cancel
hok = uicontrol('Style','pushbutton','String','OK','Position',[600,40,80,25],'Callback',{@OKFcn},'BackgroundColor',bgclr);
hcancel = uicontrol('Style','pushbutton','String','Cancel','Position',[510,40,80,25],'Callback',{@CancelFcn},'BackgroundColor',bgclr);

hfctlbl = uicontrol('Style','text','String','Cost function:','HorizontalAlignment','left','BackgroundColor',bgclr);

% Threshold text
hthresholdlbl = uicontrol('Style','text','String','Threshold:','HorizontalAlignment','left','BackgroundColor',bgclr);
hthreshold = uicontrol('Style','edit','String',num2str(threshold),'BackgroundColor','white','Callback',{@ThresholdFcn});

% Order slider
horderlbl = uicontrol('Style','text','String','Polynomial order:','HorizontalAlignment','left','BackgroundColor',bgclr);
horder = uicontrol('Style','slider','SliderStep',[0.5 0.5],'Min',0,'Max',10,'Value',order,'SliderStep',[0.1 0.1],'Callback',{@OrderFcn});
horderval = uicontrol('Style','text','String',num2str(order),'BackgroundColor',bgclr);

% Move the GUI to the center of the screen
movegui(hwin,'center');

% Plot a first estimation
[ztmp,atmp,ittmp,order,threshold,costfunction] = compute(n,y,order,threshold,costfunction);

% Make the GUI visible
set(hwin,'Visible','on');

% Callback functions

function CancelFcn(source,eventdata)
% Just close the window
uiresume(gcbf);
close(hwin);
end

function OKFcn(source,eventdata)
% Return the current estimation and close the window
z = ztmp;
a = atmp;
it = ittmp;
ord = order;
s = threshold;
fct = costfunction;
uiresume(gcbf);
close(hwin);
end

function CostFunctionFcn(source,eventdata)
% Change cost function
cf = get(hfct,'Value');
if cf == 1,
costfunction = 'sh';
elseif cf == 2,
costfunction = 'ah';
elseif cf == 3,
costfunction = 'stq';
elseif cf == 4,
costfunction = 'atq';
end
[ztmp,atmp,ittmp,ord,s,fct] = compute(n,y,order,threshold,costfunction);
end

function OrderFcn(source,eventdata)
% Change order
order = get(horder,'Value');
set(horderval,'String',num2str(order));
[ztmp,atmp,ittmp,ord,s,fct] = compute(n,y,order,threshold,costfunction);
end

function ThresholdFcn(source,eventdata)
% Change threshold
threshold = get(hthreshold,'String');
threshold = str2double(threshold);
[ztmp,atmp,ittmp,ord,s,fct] = compute(n,y,order,threshold,costfunction);
end

function [ztmp,atmp,ittmp,order,threshold,costfunction] = compute(n,y,order,threshold,costfunction)
% Compute and plot an estimation (need to sort the data)
[ztmp,atmp,ittmp,order,threshold,costfunction] = backcor(n,y,order,threshold,costfunction);
[~,i] = sort(n);
plot(n(i),y(i),'b-',n(i),ztmp(i),'r-');
end

function WinResizeFcn(source,eventdata)
% Resize the window
pos = get(hwin,'Position');
w = pos(3);
h = pos(4);
if w>400 && h>100,
set(haxes,'Position',[40,40,w-320,h-70]);
end;
set(hok,'Position',[w-90,30,80,25]);
set(hcancel,'Position',[w-180,30,80,25]);
set(hfctlbl,'Position',[w-240,h-30,220,20]);
set(hfct,'Position',[w-240,h-50,220,25]);
set(hthresholdlbl,'Position',[w-240,h-80,220,20]);
set(hthreshold,'Position',[w-240,h-100,220,20]);
set(horderlbl,'Position',[w-240,h-130,220,20]);
set(horder,'Position',[w-210,h-150,190,20]);
set(horderval,'Position',[w-240,h-150,20,20]);
end

uiwait(gcf);

end```