%		TDeious Fourier Transformation
%Number of pairs of Fourier coefficients to use in reconstruction
npairs=20
fcs=npairs+1
% Get some data
load testspec;
d=H_XV';
dp=size(d,1)
w=[1100:2:2498]';
%use this spectrum for FT
% calculates FT using equations on p186 of McClure's 
% article in Handbook of NIRA
hold on
plot(w,d,'k');
% 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;
plot(w,e1,'b');
pause
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

% invert the transform
%make unrequired coefficients 0
fc(fcs:k)=0;
fc(k+fcs:n)=0;
xfit=X*fc;
%Un-tilt the spectrum
xfits=xfit'+e';
ds=e1-xfit;
plot(w,ds,'g');
pause
plot(w,xfits,'r');
hold off
pause
figure
plot(w,ds,'m');