function [srf, s] = int2srf_nozf(isrf, s)

% function [srf, s] = int2srf_nozf(isrf, s)
%
% Convert double-sided AIRS test run interferogram to spectral SRF.
% This version does not do zero-filling.
%
% Input:
%    isrf = array (npts x nchan) of interferograms
%    s = matlab structure including of the following:
%       s.testnum = test number
%       s.module = module ID
%       s.vlaser = Interferometer HeNe laser frequency (cm-1)
%       s.dxlaser = Interferometer HeNe laser wavelength (cm)
%       s.npts = number of points in interferogram
%       s.fringes = number of HeNe laser fringes between samples
%       s.hfov = Interferometer input Half Field Of View angle (radians)
%       s.nchan = number of channels in module
%       s.chanid = channel IDs
%       s.exp_vmin = expected approx min channel freq in module
%       s.exp_vmax = expected approx max channel freq in module
%       s.dd = Interferometer sample step size (cm)
%
% Output:
%    srf = array (256 x nchan) of response functions
%    s = the input "s" structure plus the additional following items:
%       s.freq = frequency points for SRF
%       s.opd = Optical Path Difference, true (cm)
%       s.opdzf = Optical Path Difference, Zero-Filled (cm)
%       s.indok = indices of channels successfully coverted to SRFs
%       s.noisearray = measure of noise in far left & right wing of SRF
%       s.imaxarray = max value appearing in the interferogram
%       s.linearray = slope & y-intercept of straight line fit to infgrm
%       s.fmaxarray = frequency of max point in SRF; may be bogus if aliased
%       s.larray = interferogram phase correction "l"
%

% Created by Scott Hannon, 2 March 1999 (based on scripts rd_infgrm & int2srf)
% Modified 8 Mar 99, HM; split off SRF code with structure parameters
% Last updated 26 April by Scott.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%
% Unpack input structure fields needed for this routine
% Note: this is done just so that the old code below could be used
% without changing all these vars to s.* vars.
%
npts = s.npts;
nchan = s.nchan;
fringes = s.fringes;
exp_vmin = s.exp_vmin;
exp_vmax = s.exp_vmax;
testnum = s.testnum;
module = s.module;
vlaser=s.vlaser;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% old code below %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: some unnecessary lines in the old code have been deleted

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% A brief description of the phase correction technique employed %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This routine does a phase correction to the ifft of the measured
% interferogram, and seems to be similar to what JPL calls the "phs"
% spectra.  This version is based on the equation
%
%    Bm(v)= Bt(v) * exp(-i*2*pi*v*l),
%
% where "v" is wavenumber and "l" is the distance the measurement missed
% the origin during zero-crossing (modulo some unknown whole number of
% mirror movement steps).  "Bm" is the spectra from the fft of the measured
% interferogram, and unknown "Bt" is the "true" spectra (ie without phase
% error).  If there is no phase error, then
%
%    imag(Bt) = 0, and Bt = real(Bt)
%
% Therefore,
%
%    real(Bm) + i*imag(Bm) = Bt*cos(w) - i*Bt*sin(w), with w=2*pi*v*l
%
% which can be re-arranged to give
%
%    tan(w) = -imag(Bm)/real(Bm)
%
% This can be solved for "w", and thus "l".  "Bm" is (or at least should be)
% approximately zero everywhere except close to  the channel center frequency.
% This program solves for the optimal phase correction l for each channel.
%
%%%%%%%%%%%%%%% end of phase correction description %%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Convert interferogram to SRF %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Check that npts is OK
% npts
if (npts < 2^10 | npts > 2^15)
   npts
   error('npts is out of range; expecting 2^10 < npts < 2^15')
end

% Create srf and freq arrays for main output
srf=zeros(256,nchan);
freq=zeros(256,nchan);

% Laser wavelength (cm)
dxlaser=1./vlaser;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check to see if data might be aliased
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the first alias freq
alias_freq=vlaser/(fringes*2);
alias=0;
if (exp_vmax > alias_freq)
   alias=1;
   msg=['Note: will process this data assuming it is the first alias'];
   disp(msg)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Tweak min/max expected chan freqs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
chanavg=round(0.5*(exp_vmin + exp_vmax));  % Expected average freq in array
chanmin=round(exp_vmin - chanavg/170);
chanmax=round(exp_vmax + chanavg/170);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create full freq array for points (double-sided)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nhalf1=floor(0.5*npts + 0.00001);  % half, rounded down
nhalf2=npts - nhalf1;              % half, rounded up
%
totdist=(npts - 1)*dxlaser*fringes;  % total distance (max OPD is half this)
dv = 1/totdist;   % freq spacing = 1/( 2 * max opd)
% Note: the freq array ordering should be inside-out
vmin=-dv*nhalf1;
vmax=dv*(nhalf2 - 1);
fzf=[0:-dv:vmin, vmax:-dv:dv];
if (alias == 1)
   fzf=2*alias_freq - fzf;
end
%
% Determine the indices spanning the range of expected channel freqs
% If the peak of the SRF is found to be outside this range it probably
% means the measurement is no good.
indallowed=find(fzf > chanmin & fzf < chanmax);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine index offsets for noise level estimates
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The estimates are done between ~10 to 15 FWHM in the wing
noise_wing_min=10;  % Distance (FWHM) from chan center to start noise est.
noise_wing_max=15;  % Distance (FWHM) from chan center to end noise est. 
indoffmin=round(chanavg*noise_wing_min/(1200*dv));
indoffmax=round(chanavg*noise_wing_max/(1200*dv));
indnoiseoffset=indoffmin:indoffmax;
clear indoffmin indoffmax chanavg

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define optical path differences
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note: these OPDs are not directly used by this program.  They are
% included only for information purposes in the output file.
opd=0.5*(npts - 1)*dxlaser*fringes;  % True OPD
opdzf=opd;  % No zero filling

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create/declare some of the arrays that will appear in the output file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create array for raw maximum interferogram value
imaxarray=zeros(1,nchan);

% Create array for freq at ifft(isrf) max
% Note that due to low freq noise etc the peak may not always be
% at the channel center freq.
fmaxarray=zeros(1,nchan);

% Create array for zero-offset distance "l"
larray=zeros(1,nchan);

% Create array for linear (sloped straight line) fit of interferogram
linearray=zeros(2,nchan);  % two terms: slope and y-intercept

% Create array for left and right estimated RMS noise level
noisearray=zeros(2,nchan);

% Create array of indices for right (ie positive freq) half of spectra
indright=(nhalf1 + 1):(npts);

xi=(1:npts)';
indbad=[];
indok=[];
% indbadl=[];

%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over the channels
%%%%%%%%%%%%%%%%%%%%%%%%
for ichan=1:nchan

   %%%%%%%%%%%%%%%%%%%%%%
   % Check data for zeros
   %%%%%%%%%%%%%%%%%%%%%%
   % If many of the points are exactly zero, the data is no good.
   nzero=length( find(isrf(:,ichan) == 0.0) );
   if (nzero < 0.001*npts)

      %%%%%%%%%%%%%%%%%%%%%
      % Convert isrf to srf
      %%%%%%%%%%%%%%%%%%%%%
      % Find max raw interferogram value
      imaxarray(ichan)=max(isrf(:,ichan));

      % Pull out the mean from the interferogram
      % sloped straight line
      coef=polyfit(xi,isrf(:,ichan),1);
      isfmean=coef(1)*xi + coef(2);
      linearray(:,ichan)=coef';
      %
      yzf=(isrf(:,ichan) - isfmean)';

      % Do the fft (or ifft)
      % Note: need to "fftshift" to re-arrange data for use with "ifft"
      yzfifft=ifft( fftshift(yzf) );

      %%%%%%%%%%%%%%%%%%%%%%%%
      % Find the absolute peak
      %%%%%%%%%%%%%%%%%%%%%%%%
      % Note: the absolute peak is not used by this program.  It is included
      % only for information purposes in the output file.
      ymax=max( abs( real(yzfifft(indright)) ) );
      imaxright=find( abs( real(yzfifft(indright)) ) == ymax);
      imax=indright(imaxright); % Position in total array
      fmax=fzf(imax);
      fmaxarray(ichan)=fmax;

      % Check that peak freq is in expected freq range
      % if (fmax < chanmin | fmax > chanmax)
      %   msg=['Channel ' int2str(ichan) ' peak freq=' int2str(round(fmax)) ...
      %      ' is outside expected ' int2str(chanmin) ' to ' ...
      %        int2str(chanmax) ' range'];
      %   disp(msg)
      % end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Find peak in expected freq range
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      ymax=max( abs( real(yzfifft(indallowed)) ) );
      imaxallowed=find( abs( real(yzfifft(indallowed)) ) == ymax);
      imax=indallowed(imaxallowed); % Position in total array
      fmax=fzf(imax);

      if ( real(yzfifft(imax)) < 0)
         % msg=['Flipping sign of fft(isrf) for channel ' int2str(ichan)];
         % disp(msg)
         yzfifft=-yzfifft;
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Solve for phase correction "l" at the peak
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      larray(ichan)=atan( -imag( yzfifft(imax) )/real( yzfifft(imax) ) )/ ...
         (2*pi*fmax);

      % Check "l" to see if it is reasonable; if not, set to zero
      if ( abs(larray(ichan)) > 1E-3 )
         % msg=['Setting phase correction to 0 for channel ' int2str(ichan)];
         % disp(msg)
         larray(ichan)=0;
         % indbadl=[indbadl, ichan];
         Bt=real( yzfifft );
      else
         % Apply the phase correction
         Bt=real( yzfifft.*exp(i*2*pi*larray(ichan)*fzf) );
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Pull out the SRF points near the peak
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Note: the peak may have moved slightly due to phase correction
      ymax=max( Bt(indallowed) );
      imaxallowed=find( Bt(indallowed) == ymax);
      imax=indallowed(imaxallowed); % Position in total array
      % Put peak in 129th position in array of 256 points
      % Note: this is in reverse order compared to JPL
      indnear=(imax-128):(imax+127);
      indnoise1=imax - indnoiseoffset;
      indnoise2=imax + indnoiseoffset;
      if ( max(indnear) > npts)
         msg=['Unable to process data for channel ' int2str(ichan)];
         disp(msg)
         indbad=[indbad, ichan];
      else
         freq(:,ichan)=fzf(indnear)';

         if ( freq(129,ichan) < chanmin | freq(129,ichan) > chanmax )
            msg=['Freq for channel ' int2str(ichan) ...
               ' is outside expected range; assuming bad data'];
            disp(msg)
            indbad=[indbad, ichan];
         else
            Bt=Bt'/ymax;
            srf(:,ichan)=Bt(indnear);
            indok=[indok, ichan];
            noisearray(1,ichan)=sqrt(mean(Bt(indnoise1).^2));
            noisearray(2,ichan)=sqrt(mean(Bt(indnoise2).^2));
         end
      end

   else
      % Data has too many zeros; skip this channel and label it "bad"
      indbad=[indbad, ichan];
   end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% New code code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Add a few new items to the "s" output structure
%
s.freq = freq;
s.opd = opd;
s.opdzf = opdzf;
s.indok = indok;
s.noisearray =noisearray;
s.imaxarray = imaxarray;
s.linearray = linearray;
s.fmaxarray = fmaxarray;
s.larray = larray;

%%%%%%%%%%%%%%%%%%%%%%%%%%% end of routine %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
