%--------------------------------------------------------------------------
% Familiarize yourself with Matlab's Fast Fourier Transform
% This routine produces a cosine and a sine function, each with a specified
% wavelength, and allows you to see the "fast Fourier transform" (fft) of 
% these functions or the sum of the functions.  This allows you to see how
% fft's work and where Matlab puts the coefficients of these periodic
% functions
%--------------------------------------------------------------------------
clear;

%---------------------------------------------------------------
% Input variables for test curves
% A, ncycles1 = amplitude and # of cycles (which controls the wavelength)  
%               of the cosine function
% B, ncycles2 = amplitude and # of cycles (which controls the wavelength) 
%               of the sine function
% h0	      = background constant
% iwhich      =1 plot results of cos function only
%             =2 plot results of sin function only
%             =3 plot results of the combined (sin + cos) function
%--------------------------------------------------------------
A=0;	ncycles1=1;
B=0;	ncycles2=1;
h0=1;
iwhich=1;               
%----------------------------------------------------------------
% Initialize/define distance, wavenumber, and wavelength array
% nx=length of distance array
% nk=nx/2 because max. wavenumber is the inverse of twice the grid
%    interval.  This is the "Nyquist" frequency (or wavenumber)
%----------------------------------------------------------------
x=[0:5:100];   			%distance array (arbitrary units)
nx=length(x);           % length of x
dx=x(2)-x(1);           %grid spacing in x

nk=floor(nx./2);			%Indice of the Nyquist wavenumber (i.e., the
                            %maximum resolvable wavenumber)
dk=2.*pi./(max(x)+dx);		%Minimum wavenumber

if (nk < nx./2);            %Now set up the wavenumber array
  k=[[0:nk] -[nk:-1:1]].*dk;	%if nx is odd
else
  k=[[0:nk] -[nk-1:-1:1]].*dk; %if nx is even
end;

%--------------------------------------------------------
% Cosine function and fft
%--------------------------------------------------------
k1=ncycles1.*dk;		%
h1=A.*cos(k1.*x)+h0;    % h is a cos function in time
fh1=fft(h1)./nx;		% 1-D Fourier Transform. Note the fft is amplified
                        % by a factor of nx so for *display only*, divide by
                        % nx so the amplitude can be related to that of the
                        % input function.  
%--------------------------------------------------------
% Sine function and fft
%--------------------------------------------------------
k2=ncycles2.*dk;
h2=B.*sin(k2.*x)+h0;		% h is a sine function of timeMi=37.214*pi/180; Md=9.766*pi/180; 

fh2=fft(h2)./nx;		% 1-D Fourier Transform
%--------------------------------------------------------
% Combined function
%--------------------------------------------------------
h3=h1+h2-h0;			%subtract 1 h0 because its in both h1 & h2
fh3=fft(h3)./nx;

%--------------------------------------------------------
% Plots
%--------------------------------------------------------
if (iwhich==1)
  h=h1;
  fh=fh1;
elseif (iwhich==2)
  h=h2;
  fh=fh2;
else
  h=h3;
  fh=fh3;
end

figure(1); clf;
subplot(311)
plot(x,h);
xlabel('Distance (arbitrary units)');
ylabel('Height (meters)');
grid on
axis([0 max(x) -2 4]);

subplot(312)
plot([1:nx], real(fh),'o-'); hold on;
plot([1:nx], imag(fh),'r*--');
xlabel('Grid number in FFT array');
ylabel('Amplitude of Fourier coefficient (*1/nx, meters)');
axis([1 nx -1 1]);
legend('real part','imaginary part','Location','North')
grid on

subplot(313)
%plot([1:nx],2.*pi./k,'bo-');
plot(k,real(fh),'o-'); hold on;
plot(k,imag(fh),'r*-'); hold on;
legend('real part','imaginary part')
xlabel('Wavenumber, k');
ylabel('Amplitude of Fourier coefficient (*1/nx, meters)');
axis([min(k) max(k) -1 1]);
grid on;

