%Convert from wavenumber to wavelength

% Get wavenumber spectrum
load HS1;
d=(S_XV(2,1:2640))'; 
% This spectrum was measured from 9090 wn to 3812 wn
% for a total of 2640 points
dp=2640;

wnmax=9090.89;  
wngap=1.92776;
wnmin=wnmax-(dp-1)*wngap;
wnr=[wnmax:-wngap:wnmin];
figure
plot(wnr',d)
set(gca,'XDir','reverse')
% calculates FT using equations on p186 of McClure's 
% article in Handbook of NIRA

% Tilt spectrum to make ends equal
dm=dp-1;
%calculate slope
s=(d(dp)-d(1))/dm;
e=s*[0:dm]';
%adjust data by slope
e1=d-e;

y=e1;
n=size(y,1);  % number of data points
k=n/2;        % half number of data points
t=[1:n]';     % n x 1 col vector 1,2,...,n
p=[1:k];      % 1 x k row vector 1,2,...,k

wpt=(2*pi/n)*(t*p);
% n x k matrix with ij'th element 2(pi)ij/n
C=cos(wpt);   % n x k matrix of cosines
S=sin(wpt(:,1:(k-1)));   
% n x k-1 matrix of sines (omit last one
% because it would be a col of zeros anyway)
X=[ones(n,1) C S];      % put together into n x 2M matrix
div=n*[1 (1/2)*ones(1,k-1) 1 (1/2)*ones(1,k-1)]'; 
% 2k x 1 col vector of divisors 
% fc are the Fourier coefficients
% order is a0,a1,...,aM,b1,...,b(k-1)
% Fred McClure orders his differently

fc=(X'*y)./div;  % This does the summations

% now we invert the transform, not at the equally spaced 
% wavenumber points we started from, but at equally spaced 
% wavelength points

% this will involve some interpolation, so we have to smooth a
% little by omitting the very high order coefficients

% Number of pairs of Fourier coefficients to use when inverting
npairs=330;    
fcs=npairs+1;

%make unrequired coefficients 0
fc(fcs:k)=0;
fc(k+fcs:n)=0;

% set up vector with points at which to invert
w=[1100:2:2498]';

% convert this to wavenumbers
wn=(1e7)./w;

% now convert these to a scale from 1 to 2640
ti=1+(wnmax-wn)/wngap;

% now create an X matrix just like the one used for
% the FT, but with different t
ni=size(ti,1);
wpt=(2*pi/n)*(ti*p);
C=cos(wpt);   
S=sin(wpt(:,1:(k-1)));   
Xi=[ones(ni,1) C S];   

% now invert, but using Xi instead of X
xfit=Xi*fc;

figure
plot(w,xfit,'g');
pause;
hold on

% get dispersive spectrum and plot for comparison
load HS2;
dd=(S_XV(2,:))'; 
plot(w,dd,'r');
pause
% Un-tilt the spectrum (Not really necessary)
% be careful, need to convert to 1-2640 scale on which
% linear baseline was subtracted
ei=s*(ti-1);
xfits=xfit'+ei';
plot(w,xfits);

hold off