resample_rev2(sigt, masst, sigma0, skew0, pfit, minHW, Gnoise)

function [sigr, massr, rrr,ii, nsigma, nskew] = resample_rev2(sigt, masst, sigma0, skew0, pfit, minHW, Gnoise)
%   
% Function resamples TOF and m/z axis and integrates signal with the step
% determined by the ratio of predicted HWFM for the peak outside
% mass focusing range (given by quadratic polynomial fit) to the
% intitial FWHM (left half) of the peak in the mass focusing range
% as measured by "peakwidth_pf"
%
% INPUT:
% 
% sigt - raw equisampled TOF signal;
% masst - mass axis, corresponding to sigt;
% sigma0 - value of FWHM for the instrumental function peak in the
%          mass focusing range (in time points, can be fractional);
% skew0 - right peak half skew (can be negative)
% pfit - is the array of polyfit coefficients [p1,p2,p3] for the
% expected peak HWFM  outside mass focusing range;
% minHW - sufficient minimum peak point density (4.5 by default)
% Gnoise - set to 1 in case of Gaussian noise (e.g., 0 in sace of
% derivative noise)
%
% Output:
%
% sigr - resampled, integrated TOF signal;
% massr - resampled mass axis, corresponding to sigr
% rrr - resampling rate curve;
% ii - resampled time;
% nsigma - new sigma, if constant resampling was used;
% nskew - new skew, if constant resampling was used;
%
% USAGE:
% [sigr, massr, rrr,ii, nsigma, nskew] = resample_rev2(sigt, masst, sigma0, skew0, pfit, minHW, Gnoise);
% Dependency: none
%

if (nargin < 2) | (nargin < 1)
  error('Requires both the data array and the peak position as an input');
  sigma = NaN;
end

nsigma=sigma0; nskew=skew0;

t=1:length(sigt);
t=t'; tb=t;
rr=(pfit(1)*t.*t+pfit(2)*t+pfit(3))/sigma0;
rr0=rr;
stp=round(sigma0/minHW);
if stp >= 3
  % 'doing_initial_constant_resampling'
  if (stp/2-fix(stp/2))==0 stp=stp-1; end;
  trs=1:stp:length(t);
  k=1;sigr_stp=zeros(length(trs),1);
  ks=(stp-1)/2;
  sigr_stp(1)=sum(sigt(1:ks));
  sigr_stp(end)=sum(sigt(trs(end)-ks:trs(end)));
  
  for jk=trs(2):stp:trs(end-1)
    k=k+1; 
    sigr_stp(k)= sum(sigt(jk-ks:jk+ks));
  end;
  stp;  % constant IDS step
  nsigma=sigma0/stp; nskew=skew0/stp;
  masst=masst(trs);
  sigt=sigr_stp;
  rr0=rr0(trs);
  rr=rr(trs);
  t=t(trs);
end;

t1=(1:length(t))'; % assuming that length(sigt)=length(t)
[minrr0, iminrr0]=min(rr0);
%plot(rr);
il=find((rr-floor(rr))<0.5);
rr(il)=floor(rr(il));

ih=setdiff(t1,il);

rr(ih)=ceil(rr(ih));

[minrr, iminrr]=min(rr);

ir=find(rr>=2); ir=ir(find(ir>iminrr));
%length(ir)

if(length(ir)<1) %% if any quadratic resampling should be done
  sigr=sigt;
  massr=masst;
  rrr=ones(length(t),1);
  ii=t;
  
else
  
  % 'onset'
  ir(1)*stp  % onset of sampling in "original" time
  %'max_rate'
  mr=max(rr);  % maximum sampling rate for this record
  %plot(rr)
  
  sigr=sigt(1:ir(1)-2); sigr(end)=sigr(end)+0.5*sigt(ir(1)-1);
  %%%!!! need to weigh original points between 1 and 1.5
  massr=masst(1:ir(1)-2);
  rrr=rr(1:ir(1)-2);
  
  ii=[(1:(ir(1)-2))'];
  %'last_ii'
  ii(end);
  rsp=ir(1);
  rst=ir(1);
  
  for jj=(rr(ir(1))+1):mr
    %  jj
    ich=find(rr>=jj); ich=ich(find(ich>ir(1)));
    rsp=[rsp,ich(1)];
    % 'first_ii'
    % correct for integer sampling
    rr((ii(end)+jj-1):rsp(end-1))=rr(rsp(end-1));
    %ii(end)+jj-1
    rst=[rst, ii(end)+jj-1];
    ii=[ii;t1((ii(end)+jj-1):(jj-1):(rsp(end)-jj))];
    % 'last_ii'
    % ii(end)
  end;
  
  %'first_ii'
  % back correct for integer jump lagging
  rr((ii(end)+mr):rsp(end))=rr(rsp(end));
  %ii(end)+mr
  rst=[rst,ii(end)+mr]; rst=rst(2:end);
  
  ii=[ii;t1((ii(end)+mr):mr:(t1(end)-ceil(mr/2)))];
  
  %'last_ii'
  %ii(end);
  %%rsp  % rounded resampling start for different rates
  %%rst  % start resampling corrected for integer jumps
  %% [rr(rst), rr(rst-1), rr(rst+1)] % check the edge points
  %plot(ii, 'k.');
  
  irs=ii(rst(1):end); % resampling time index
  massr=[massr;masst(irs)]; % resampled mass
  rrr=[rrr;rr(irs);rr(irs(end))]; % resampling rate
  
  for ri=1:(length(irs)-1)
    rri=rr(irs(ri));
    if (rri/2-fix(rri/2))==0 
      k=rri/2;
      sri=sum(sigt(irs(ri)-k+1:irs(ri)+k-1))+0.5*sigt(irs(ri)-k);
      rri1=rr(irs(ri+1));
      if (rri1/2-fix(rri1/2))~=0 
	sri=sri+sigt(irs(ri)+k); 
      else sri=sri+0.5*sigt(irs(ri)+k);
      end;
    else
      k=(rri-1)/2;
      sri=sum(sigt(irs(ri)-k:irs(ri)+k));
      rri1=rr(irs(ri+1));
      if (rri1/2-fix(rri1/2))==0 
	sri=sri+0.5*sigt(irs(ri)+k+1); 
      end;
    end;
    
    sr(ri)=sri;  
  end;   % loop through resampled signal
  sigr=[sigr;sr'];  %% main work is done for signal intensity integration
  %% resample last point assuming that next rri1 is the same
  ri=length(irs); rri=rr(irs(ri));
  if (rri/2-fix(rri/2))==0 
    k=rri/2;
    sri=sum(sigt(irs(ri)-k+1:irs(ri)+k-1));
    sri=sri+0.5*sigt(irs(ri)-k)+0.5*sigt(irs(ri)+k);
  else
    k=(rri-1)/2;
    sri=sum(sigt(irs(ri)-k:irs(ri)+k));
  end;
  sigr=[sigr;sri;mean(sigr(end-k:end))]; %% add last resampled point
  %% assumes that last point is a mean of (end-k): 
  %% one point lost
  
  %%% add the last time point of the data as "mean" of what is left
  ii=[ii;t1(end)]; 
  rrr=[rrr;rrr(end)];
  massr=[massr; masst(ii(end-1));masst(t1(end))];
  sigr=[sigr; mean(sigt(ii(end-1):ii(end)))];
  
  %%% scale by sqrt(rate) for Gaussian noise
  if Gnoise==1    
  sigr((ir(1)-2):end)=sigr((ir(1)-2):end)./((rrr((ir(1)-2):end)).^(0.5));    
  end;
  
  %%%%% try to correct for ramps %%%%%%
  rmpi=find(diff(rrr)==1); rmpi=rmpi+1; % ramp positions
  rmpa=sigr(rmpi)-sigr(rmpi-2); % ramp amplitudes
  %determine noise amplitude
  [y,mi]=min(sigr(rmpi(1):(end-2*mr-1)));
  pst=rmpi(1)+mi;
  nsa=2*std(sigr((pst-2*rrr(pst)):(pst+2*rrr(pst))));
  i1=find(rmpa>nsa & ((rmpa+nsa)>(sigr(rmpi+1)-sigr(rmpi-1)))); 
  ii(rmpi(i1)); % ramps to correct
  for cr=1:length(i1)
    si=rmpi(i1(cr));  % cor-ramp positions in resampled signal
    trmp=ii(rmpi(i1(cr))); % cor-ramp positions in original time
    trmp1=ii(rmpi(i1(cr))-2);
    %iminrr0
    rrmp=rrr(rmpi(i1(cr))); % resampling rate at the ramp
    lcr=fix(length(find(rr0(trmp:end)<rrmp))/rrmp); 
    % length to the right
    %rr0(trmp+lcr*rrmp)
    %rrr(si+lcr)
    if (lcr<1) lcr=1; end; %% avoid 0-division 11/07/07
    wcr=1./lcr; % weight step to the right
    lcr=floor(lcr-2*nsa/(wcr*rmpa(i1(cr)))); 
    wcl=1; % initialize
    % limited length to the right
    if (lcr>0) wcr=1./abs(lcr);  end; % if negative??
    % rr0(trmp+abs(lcr)*rrmp)
    %rrr(si+abs(lcr))
    lcl=fix(length(find(rr0(iminrr0:trmp1)>(rrmp-1)))/(rrmp-1));
    %rr0(trmp1-lcl*(rrmp-1))
    %rrr(si-2-lcl)
    if (lcl<1) lcl=1; end; %% avoid 0-division 11/07/07
    if (lcl>0)  % if negative??
    wcl=1./lcl; % weight step to the left 
    lcl=floor(lcl-2*nsa/(wcl*rmpa(i1(cr)))); 
    % limited length to the left
    wcl=1./abs(lcl);
    end;
    %rr0(trmp1-abs(lcl)*(rrmp-1))
    %rrr(si-2-abs(lcl))
    rrl=length(find(rrr==rrmp));
    rll=length(find(rrr==(rrmp-1)));
    
    if (abs(lcr)>0.5*rrl) lcr=floor(0.5*rrl); wcr=1./lcr; end;
    
    dcr=0.5*rmpa(i1(cr))*(1-wcr*(1:abs(lcr))');
    sigr(si:(si+abs(lcr)-1))=sigr(si:(si+abs(lcr)-1))-dcr;
    
    if (abs(lcl)>0.5*rll) lcl=floor(rll/2); wcl=1./lcl; end;
    
    dcl=0.5*rmpa(i1(cr))*(1-wcl*(1:abs(lcl))');
    sigr((si-2):-1:(si-2-abs(lcl)+1))=sigr((si-2):-1:(si-2-abs(lcl)+1))+dcl;
    
  end; % loop through corrected ramps
  ii=t(ii); 
  %ir(1)
  %ii(ir(1))
  ii(ir(1):(end-1))=ii((ir(1)+1):end); %% correct for 1-point resampling onset shift 11/06/07
  ii(end)=tb(end); %% correct for copy of end-value 11/19/07
end; % if resampling length is nonzero (safety check)

return;