

% Whittaker Smoother method for GNAT data

% Author: Marshall Smith

% References: 1) Paul H. C. Eilers Anal. Chem. 2003, 75, 14, 3631–3636
%             2) J. Carlos Cobas, Michael A. Bernstein, Manuel Martin-Pastor, Pablo Garcia Tahoces J. Mag. Reson. 2006, 183, 145-151


function whsm(~,~)
%% INPUT values for the user

prompt={'enter required smoothings factor'};
title={'Whittaker Smoother Parameters'};
fieldsize = [1 45];

l={'10'}; % Smoothing factor
lw=1;
answer = inputdlg('Enter Smoothing factor',...
             'Whittaker Smoother', [1],l);
lambda=str2double(answer{1});

doplot=0; % set to 1 to show comparison of original and corrected spectrum

% Get the data from the main GUI
GNATD.data=findobj(allchild(0), 'tag', 'MainFigure');
GNATD.NmrData=guidata(GNATD.data(1)); % Get the NmrData structure from the main window
spectra=real(GNATD.NmrData.SPECTRA);
[m,n]=size(spectra);
cSpectrum=zeros(m,n); %preallocate the array
psw=GNATD.NmrData.sw*GNATD.NmrData.sfrq/GNATD.NmrData.np;
lwp=round(lw/psw/2); 

for k=1:n
spectrum=spectra(:,k);

if GNATD.NmrData.th==0
    errordlg('Please set a threshold')
end 

 peak=peakpick_mn(spectrum,GNATD.NmrData.th);  % take user input for threshold
 peak=table2array(struct2table(peak));

 mask=zeros(1,length(spectrum));
 for i=1:length(peak) 
     mask(peak(i,2)-lwp:peak(i,3)+lwp)=1; 
 end
  
%% estimate the baseline
if(size(spectrum,1) == 1)
  spectrum = spectrum';
end
if(size(mask,1) == 1)
  mask = mask';
end

%invert logical value of mask
maskInv = sparse(ones(size(mask)) - mask);

% create the D vector, the finite difference operator
N = length(spectrum);
D = sparse(diag(ones(N, 1), 0) - diag(ones(N-1, 1), -1));

% create the W matrix using sparse diagonal
W = spdiags(maskInv,0,N,N);

% create the D' * D matrix
DtD = D' * D;

% set up the linear system to solve
A = W + lambda * DtD;
b = W*spectrum;


% solve the system efficiently using iterative solver (e.g., pcg)
 tolerance = 1e-6;
 maxIterations = 1e6;
 base3 = bicgstab(A, b, tolerance, maxIterations); 
cSpectrum(:,k) = spectrum - base3;

end 

if doplot==1
figure();
  hold on;
  plot(real(spectrum));
  scaledMask = mask* max(abs((spectrum)));
  plot(real(scaledMask),'g-')
  plot(real(base3),'r-')
plot(cSpectrum)
  hold off
end 

%% UPDATE THE MAIN GUI DATA 
% DOES NOT UPDATE PLOT AUTOMATICALLY PRESS SPECTRUM/Full UNDER PLOT IN GNAT TO 
% UPDATE THE SPECTRUM
 
 GNATD.NmrData.SPECTRA=cSpectrum;
 NmrData=GNATD.NmrData;

 guidata(GNATD.data(1),NmrData); %Updates the data in Main figure but doesn't plot press full to update figure
end

    