function [MZOUT,YOUT] = msresample(MZ,Y,N,varargin)
%MSRESAMPLE resamples a mass spectrometry signal
%
%   [MZOUT,YOUT] = MSRESAMPLE(MZ,Y,N) resamples a mass spectrometry signal.
%   Y are the ion intensities and MZ is the mass-charge vector. The output
%   spectrum will have N samples where the spacing increases linearly
%   within the range [min(MZ) max(MZ)]. I.e., MZ is assumed to be a
%   quadratic function of its index. When input arguments are set such that
%   down-sampling takes place, MSRESAMPLE first applies a lowpass filter to
%   minimize aliasing.
%
%   For the antialias filter a linear-phase FIR filter is designed using
%   least-squares error minimization. The cutoff frequency is set by the
%   largest down-sampling ratio when comparing the same regions of the MZ
%   and MZOUT vectors.
%
%   Y can be a matrix with several spectrograms, all sharing the same MZ
%   scale. 
%
%   MSRESAMPLE(...,'UNIFORM',U) forces the MZ vector to be uniformly
%   spaced when U is TRUE. U defaults to FALSE.
%
%   MSRESAMPLE(...,'RANGE',R) uses a 1-by-2 vector with the mass-charge
%   range for the desired output signal. R values must be within [min(MZ)
%   max(MZ)]. R defaults to [min(MZ) max(MZ)].  
%
%   MSRESAMPLE(...,'MISSING',M) analyzes an MZ vector for dropped samples
%   when M is TRUE. M defaults to FALSE. Note: checking for dropped samples
%   can be time consuming and might not be worthwhile if the down-sample
%   factor is large. Dropped samples can only be recovered if the original
%   MZ values follow a linear or a quadratic function of the sample index.
%
%   MSRESAMPLE(...,'WINDOW',W) sets the window used during the filter
%   design process. W defaults to 'Flattop'. Other options are 'Blackman',
%   'Hamming', and 'Hanning'.
%
%   MSRESAMPLE(...,'CUTOFF',F) sets the cutoff frequency. F can be a scalar
%   value between 0 and 1. 1 corresponds to the Nyquist frequency or half
%   the sampling frequency. By default, MSRESAMPLE automatically estimates 
%   F by inspecting MZ and MZOUT; however, F can be underestimated if MZ
%   presents anomalies.
%
%   MSRESAMPLE(...,'SHOWPLOT',SP) plots the original and the resampled
%   spectra. When SP is TRUE, the first spectrogram in Y is used. If
%   MSRESAMPLE is called without output arguments, a plot will be shown
%   unless SP is false. SP may also contain an index to one of the
%   spectrograms in Y.
%  
%   Example:
%
%     load sample_hi_res
%     MZ = MZ_hi_res;
%     Y = Y_hi_res;
%    
%     % Resample the spectrogram by a 1/2 factor.
%     N = numel(MZ)/2  % find out the number of required samples
%     [mz1,y1] = msresample(MZ,Y,N);
%
%     % Resample the spectrogram to have 10000 samples between 2000 and
%     % 11000 Da. and plot the original and the resampled signals.
%     R = [2000 11000]; % set the range
%     [mz2,y2] = msresample(MZ,Y,10000,'RANGE',R,'SHOWPLOT',true);
%
%   See also MSALIGN, MSBACKADJ, MSHEATMAP, MSLOWESS, MSNORM, MSSGOLAY,
%   MSVIEWER. 

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

%set defaults
checkMissing = false;
evenlySpaced = false;
windowType = 'flattop';
cutoffGiven = false;

if nargout == 0
    plotId = 1; 
else
    plotId = 0;
end

% validate required inputs 
if size(MZ,1)~=size(Y,1)
    error('Bioinfo:resample:differentSize','MZ and Y must be the same size.')
end
if ~isnumeric(Y) || ~isreal(Y) 
   error('Bioinfo:resample:IntensityNotNumericAndReal','Ion intensities must be numeric and real.') 
end
if ~isnumeric(MZ) || ~isreal(MZ) || ~isvector(MZ)
   error('Bioinfo:resample:MZNotNumericAndReal','MZ scale must be a numeric and real vector.') 
end

R = [max(0,min(MZ)) max(MZ)];

if ~isscalar(N)
     error('Bioinfo:resample:badNValue','Number of samples (N) must be a scalar.')
end

if  nargin > 3
    if rem(nargin,2) == 0
        error('Bioinfo:resample:IncorrectNumberOfArguments',...
            'Incorrect number of arguments to %s.',mfilename);
    end
    okargs = {'range','uniform','missing','window','showplot','cutoff'};
    for j=1:2:nargin-3
        pname = varargin{j};
        pval = varargin{j+1};
        k = strmatch(lower(pname), okargs);%#ok
        if isempty(k)
            error('Bioinfo:resample:UnknownParameterName',...
                'Unknown parameter name: %s.',pname);
        elseif length(k)>1
            error('Bioinfo:resample:AmbiguousParameterName',...
                'Ambiguous parameter name: %s.',pname);
        else
            switch(k)
                case 1 % range
                    R = pval;
                    if numel(R)~=2 || diff(R)<=0 || R(1)<min(MZ) || R(2)>max(MZ)
                        error('Bioinfo:resample:badRange',...
                            'The range R must be provided in the form [MZmin, MZmax]')
                    end
                    R(1) = max(R(1),0);
                case 2 % uniform
                    evenlySpaced = opttf(pval);
                    if isempty(evenlySpaced)
                        error('Bioinfo:resample:InputOptionNotLogical',...
                            'UNIFORM must be a logical value, true or false.')
                    end
                case 3 % missing
                    checkMissing = opttf(pval);
                    if isempty(checkMissing)
                        error('Bioinfo:resample:InputOptionNotLogical',...
                            'MISSING must be a logical value, true or false.')
                    end
                case 4 % window
                    windowTypes = {'flattop','blackman','hamming','hanning','kaiser'};
                    windowType = strmatch(lower(pval),windowTypes); %#ok
                    if isempty(windowType)
                        error('Bioinfo:resample:NotValidWindow',...
                              'Not a valid window for filter design.')
                    end
                 case 5 % show
                    if opttf(pval) 
                        if isnumeric(pval)
                            if isscalar(pval)
                                plotId = double(pval); 
                            else
                                plotId = 1;
                                warning('Bioinfo:msresample: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        
                 case 6 % cutoff
                    F = pval(1);
                    cutoffGiven = true;
                    if F<=0 || F>1
                        error('Bioinfo:resample:InvalidCutoffFrequency',...
                              'Cutoff frequency must be between 0 and 1.')
                    end
            end
        end
    end
end

% make sure there are no dropped samples 
if checkMissing
    [MZ,Y] = msdroppedsamples(MZ,Y);
end

% figure out the new MZ vector
if evenlySpaced
    MZOUT = R(1)+(0:(N-1))'*diff(R)/N;
else
    MZOUT = (sqrt(R(1))+(0:(N-1))'*diff(sqrt(R))/N).^2;
end

if cutoffGiven
    DF = F;
else
    % estimate variable sampling factor (this loop is slighty faster than
    % histc)
    df = zeros(1,N);
    i = 1; j = 1;
    while j <= N
        count=0;
        while MZ(i)<MZOUT(j)
            i = i + 1;
            count = count+1;
        end
        df(j) = count;
        j = j+1;
    end
    df(1) = df(2); % the first buck contains all samples to the left of the 
                   % range 
    
    % 'df' now contains the observed variable sampling ratio at every MZ
    % value, to minimize the effect of outliers two steps are taken:
    %    1) eliminate low MZ (MZmin = 50) values because df can easily be
    %    overestimated due to the small increments when MZ follows a
    %    quadratic law. 
    %    2) fit df to a low order polynomial (polOrder = 10), so we can
    %    figure out the maximum downsampling rate without taking a likely
    %    outlier. 
    
    MZmin = 50;
    polOrder = 10;
    startAt = floor(min(N/5,find(MZOUT>MZmin,1,'first'))); 
    DF = 1/max(polyval(polyfit((startAt:N)/N,df(startAt:N),polOrder),(startAt:N)/N));
    
    % this plot is only for debugging purposes (not documented)
    plotFlag = false;
    if plotFlag 
        figure
        plot(MZOUT,df,'r'); 
        hold on
        plot(MZOUT,polyval(polyfit((startAt:N)/N,df(startAt:N),polOrder),(1:N)/N))
        title('Variable downsampling ratio')
        xlabel('Mass/Charge (M/Z)')
        hold off
    end
end

if DF >= 1  % it is UPSAMPLE, then just interpolate
    YOUT = interp1(MZ,Y,MZOUT,'spline');  
else % DF < 1  it is DOWNSAMPLE, do antialias filter before interpolate
    filtLength = max(21,round(1/DF)*4+1);
    YF = filter(msfir(filtLength,DF,windowType),1,[Y;zeros(filtLength-1,size(Y,2))]);
    YF = YF(1+(filtLength-1)/2:end-(filtLength-1)/2,:);
    YOUT = interp1(MZ,YF,MZOUT,'linear');
end

numSpectrograms = size(Y,2);
if (plotId~=0) && ~any((1:numSpectrograms)==plotId)
    warning('Bioinfo:msresample:InvalidPlotIndex',...
    'SHOWPLOT is not a valid spectrogram index, no plot was generated.')
elseif plotId > 0
    figure
    plot(MZ,Y(:,plotId),'.')
    hold on
    plot(MZOUT,YOUT(:,plotId),'r.-')
    title(sprintf('Spectrogram ID: %d  Cutoff Freq: %f ',plotId,DF))
    xlabel('Mass/Charge (M/Z)')
    ylabel('Relative Intensity')
    legend('Original samples','Up/down-sampled spectrogram')
    axis([min(MZ) max(MZ) min(Y(:,plotId)) max(Y(:,plotId))])
    grid on
    hold off
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = msfir(L,F,windowType)
% Designs a low pass anti-aliasing filter. For more advance filter design
% options see the Signal Processing Toolbox.
% Length (L) is the size of filter.
% Cut-off frequency (F) is a value between [0 1] using normalized
% frequencies. 

% Reference for window coefficients: Bruel & Kjaer, Windows to FFT Analysis
%                                   (Part I), Technical Review, No. 3, 1987 

N = L-1; % filter order

% f = firls(N,[0;F;F;1],[1;1;0;0]);
if rem(L,2)==0 % Type II linear phase filter (even)
    x = pi/2 * F * (1:2:N);
    f = F * sin(x)./x;
    f = [fliplr(f) f];
else % Type I linear phase filter (odd)
    x = pi * F * (1:N/2);
    f = F * sin(x)./x;    
    f = [fliplr(f) F f];
end

x = (2*pi/N)*(0:N);
switch windowType
    case 'hannning' % Hanning window
        w = 0.5 - 0.5*cos(x);
    case 'hamming'  % Hamming window
        w = 0.54 - 0.46*cos(x);
    case 'blackman' % Blackman window 
        w = 0.42 - 0.5*cos(x) + 0.08*cos(2*x);
    case 'flattop' % Flattop window
        w = 0.2156 - 0.4160*cos(x) + 0.2781*cos(2*x) - 0.0836*cos(3*x) +...
            0.0069*cos(4*x);
    case 'kaiser'  % Kaiser window (undocumented)
        w = kaiser(N+1,7.8562)'; % only valid when SPT is present 
end

f = f.* w; % Apply windowing
f = f/sum(f);