function Y = mssgolay(MZ,Y,varargin)
%MSSGOLAY provides least-squares polynomial smoothing of mass spectrometry
%signals 
%
%   YOUT = MSSGOLAY(MZ,Y) smoothes the mass spectrometry signal Y using the
%   least-squares digital polynomial filters (Savitzky and Golay filters).
%   The default span (or frame) is 15 samples. MZ is the mass-charge column
%   vector with the same size as the vector of ion intensity samples Y (or
%   spectrum). Y can be a matrix with several spectrograms, all sharing the
%   same MZ scale.
%
%   MSSGOLAY(...,'SPAN',S) allows you to modify the frame size for the
%   smoothing function. If S is greater than 1, then the window is of size
%   SP (in samples, i.e., independently of the MZ vector). Higher values
%   will smooth the signal more, at the expense of computation time. If S
%   is less than 1, then the window size is taken to be a fraction of the
%   number of points in the data; e.g., for S = 0.05, the window size equal
%   to 5% of the number of points in MZ. 
%
%   Notes:
%   1) The original Savitzky and Golay algorithm assumes a uniformly
%   spaced MZ vector; MSSGOLAY also allows one that is not uniformly
%   spaced. Therefore, the sliding frame for smoothing is centered using
%   the closest samples in terms of the MZ value and not in terms of the MZ
%   index.
%   2) When the MZ vector does not have repeated values or NaNs, the
%   algorithm will be approximately two times faster.
%   3) When the MZ vector is evenly spaced, the least-squares fitting is
%   performed once and the spectrum is filtered with the same coefficients,
%   speeding up the algorithm considerably.
%   4) If the MZ vector is evenly spaced and S is even, S is automatically
%   incremented by 1 to include both edge samples in the current frame. 
%
%   MSSGOLAY(...,'DEGREE',D) sets the D of the polynomial to be fitted to
%   all the points in the moving frame. The default DEGREE is 2. D must be
%   smaller than S. 
%
%   MSSGOLAY(...,'SHOWPLOT',SP) plots the smoothed spectra over the
%   original. When SP is TRUE, the first spectrogram in Y is used. The
%   default is FALSE. No plot is generated unless MSSGOLAY is called
%   without output arguments. SP can also contain an index to one of the
%   spectrograms in Y.  
%
%   Example: 
%
%       load sample_lo_res
%
%       % Smooth a group of spectrograms.
%       YS = mssgolay(MZ_lo_res,Y_lo_res);
%
%       % Plot the third spectrogram in Y_lo_res and its smoothed signal.
%       mssgolay(MZ_lo_res,Y_lo_res,'SHOWPLOT',3);
%
%   See also MSALIGN, MSBACKADJ, MSHEATMAP, MSLOWESS, MSNORM, MSRESAMPLE,
%   MSVIEWER.

% Copyright 2003-2004 The MathWorks, Inc.
% $Revision: 1.1.12.1 $   $Date: 2004/12/24 20:43:30 $


% set default parameters
span = 15;
degree = 2; 
if nargout == 1
    plotId = 0; 
else
    plotId = 1;
end

% deal with the various inputs
if nargin > 2
    if rem(nargin,2) == 1
        error('Bioinfo:mssgolay:IncorrectNumberOfArguments',...
            'Incorrect number of arguments to %s.',mfilename);
    end
    okargs = {'span','degree','showplot'};
    for j=1:2:nargin-2
        pname = varargin{j};
        pval = varargin{j+1};
        k = strmatch(lower(pname), okargs);%#ok
        if isempty(k)
            error('Bioinfo:mssgolay:UnknownParameterName',...
                'Unknown parameter name: %s.',pname);
        elseif length(k)>1
            error('Bioinfo:mssgolay:AmbiguousParameterName',...
                'Ambiguous parameter name: %s.',pname);
        else
            switch(k)
                case 1  % span
                    if ~isnumeric(pval) || ~isscalar(pval)
                        error('Bioinfo:mssgolay:LowessSpanNumeric',...
                            'SPAN must be a numeric value.');
                    end
                    span = pval;
                    if span <= 0
                        error('Bioinfo:mssgolay:LowessSpanPositive',...
                            'SPAN must be a positive value.');
                    end
                    if span >= numel(MZ)
                        error('Bioinfo:mssgolay:LowessSpanLessThanX',...
                              'SPAN must be less than the number of elements in Y.');
                    end
                case 2 % degree
                    if ~isnumeric(pval) || ~isscalar(pval)
                        error('Bioinfo:mssgolay:DegreeNumeric',...
                              'DEGREE must be a numeric value.');
                    end
                    degree = pval;
                    if degree<0
                        error('Bioinfo:mssgolay:LowessBadNumericDegree',...
                              'DEGREE should be a positive integer.');
                    end
                case 3 % show
                    if opttf(pval) 
                        if isnumeric(pval)
                            if isscalar(pval)
                                plotId = double(pval); 
                            else
                                plotId = 1;
                                warning('Bioinfo:mssgolay:SPNoScalar',...
                                'SHOWPLOT must be a single index to one of the spectrograms in Y or\na logical. Plotting the first column in Y only.')
                            end 
                        else
                            plotId = 1;
                        end
                    else
                        plotId = 0;
                    end
            end
        end
    end
end


% validate MZ and Y
if size(MZ,2) ~= 1
   error('Bioinfo:mssgolay:OnlyOneMZScale','MZ scale must be a single column.')
end
if ~isnumeric(Y) || ~isreal(Y)
   error('Bioinfo:mssgolay:IntensityNotNumericAndReal',...
         'Ion intensities must be numeric and real.') 
end
if ~isnumeric(MZ) || ~isreal(MZ)
   error('Bioinfo:mssgolay:MZNotNumericAndReal','MZ scale must be numeric and real.') 
end
if size(MZ,1) ~= size(Y,1)
   error('Bioinfo:mssgolay:NotEqualNumberOfSamples',...
         ['The ion intensity vector(s) must have the\n'...
          'same number of samples as the MZ scale.'])
end

numSpectrograms = size(Y,2);

if (plotId~=0) && ~any((1:numSpectrograms)==plotId)
    warning('Bioinfo:mssgolay:InvalidPlotIndex',...
    'SHOWPLOT is not a valid spectrogram index, no plot was generated.')
end

% iterate for every spectrogram
for ns = 1:numSpectrograms
if nargout>0 || (ns == plotId)
    
    if (ns == plotId)
        yo = Y(:,ns);
    end
    
    Y(:,ns) = masmooth(MZ,Y(:,ns),span,'sgolay',degree);
        
    if (ns == plotId)
        figure
        plot(MZ,[yo Y(:,ns)])
        title(sprintf('Spectrogram ID: %d',ns));
        xlabel('Mass/Charge (M/Z)')
        ylabel('Relative Intensity')
        legend('Original spectrogram','Smoothed spectrogram')
        axis([min(MZ) max(MZ) min(Y(:,ns)) max(Y(:,ns))])
        grid on
        hold off
    end
    
end % if nargout>0 || (ns == plotId)    
end % for ns = 1:numSpectrograms