%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  FFT demos:  C. Johnson for 450, Oct 2014
%
%  DEMONSTRATES the single approach that works for FFTs of either even or odd # of points
%  http://blogs.uoregon.edu/seis/wiki/unpacking-the-matlab-fft/
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% clean up
clear;
close all;

% Edit nx to be 256 or 257 to demonstrate that this works for either even
% or odd # points

nx=257;             % number of points; for a time series call this nt
dx=2;               % sampling interval in meters in x-direction; for a time series this is dt
lambda=64;          % wavelength; for a time series this is period
x=0:dx:dx*(nx-1);   % distances in meters; for a time seris this is time in e.g. secs
A0=1;               % amplitude
y1=A0*cos(2*pi*x/lambda);    % cosine function

%% get the fft
y1fft=fft(y1);

%%  THIS IS HOW TO ASSIGN WAVENUMBERS:  ODD or EVEN # points in your time series (2-sided spectrum)
knyq=1./(2*dx);         % Nyquist wavenumber (or freq, f if time series) = 0.5*sampling freq = 0.5*(1/dx) = 1/(2*dx); 
dk=1./(nx*dx);
kk=(0:(nx-1))*dk;       % for from k=0 to the sampling freq in steps of dk
kk(kk>knyq)=kk(kk>knyq)-knyq*2;     % now replace freqs above the nyquist with the negative freqs mirroring positive ones

%% now get 1-sided spectrum
% Manar comment: one-sided spectrum should also include the "zeroth"
% wavenumber. If you only take the positive frequencies, to properly scale
% the 1-sided spectrum you should multiply the whole thing by 2 not only
% (2:end)

kk1=kk(kk>=0);
amp1side=abs(y1fft(kk>=0));

% kk1=kk(kk>0);
% amp1side=abs(y1fft(kk>0));

%% now scale the amplitude spectrum properly
% Manar edit: the whole 1-sided spectrum (except for the zero frequency)
% should be multiplied by 2 not just half of it. nx is the length of the
% 2-sided spectrum so nx/2 is the length of the 1-sided spectrum. The
% reason to adding ceil(nx/2) is to deal with the case when there is an odd
% number of points

amp1side=amp1side/nx;
amp1side(2:ceil(nx/2))=2*amp1side(2:ceil(nx/2));  
% amp1side(2:ceil(end/2))=2*amp1side(2:ceil(end/2));  


%% make a plot and check you got the right answer
figure(99)
    subplot(211), plot(kk1,amp1side,'k-*'); 
        set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');
        xlim([0 0.05]); grid on
    subplot(212), plot(kk1,amp1side,'k-*'); 
        set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');

        
%  CHECK THAT MY PEAK IS AT THE CORRECT WAVENUMBER
fprintf('-----------------------------------------------------------\n');
fprintf('My cosine has a wavelength, lambda = 64 so k=1/64=0.0156\n\n');
ik=find(amp1side>0.1);                     % looking for the wavenumber with biggest amplitude
fprintf('Peak amplitude is = %.2f at wavenumber=%.4f \n', amp1side(ik),kk1(ik));
    