%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  FFT demos:  C. Johnson for 450, Oct 2014
%
%  DEMONSTRATES
%  1) fft in matlab for tiem series with even # of pts
%  2) what fftshift does for even and odd number of points
%  3) 2-sided spectrum, normalizing by N and how to generate wave# vector
%  4) 1-sided spectrum, normalize by 2/N and how to generate wave# vector
%  5) 2,4,5 above but for a time series with an odd # of points:  NOT SURE
%  THIS IS WORKING PROPERLY SO IGNORE FOR NOW - POSTED VERSION HAS ONLY
%  EVEN # POINTS
%
%  cells separated by "pause": hit return in command window to continue
%
%  I recommend just always making your time series have an even # of points
%  then generate the wavenumber vector using recipe in 3 or 4 above.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear;
close all;

%% DEMO 1:  the raw fft output for a cosine function

nx=256;             % 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
y2=A0*sin(2*pi*x/lambda);    % cosine function

% plot my function in the time domain
figure(1)
plot(x,y1,'k-*')
hold on
plot(x,y2,'r-*')
set(gca, 'FontSize', 18);
title('y=cos(2*pi*x/\lambda)');

pause; 

% take an fft and plot it
y1fft=fft(y1);
y2fft=fft(y2);

figure(2)
subplot(221), plot(real(y1fft)); hold on; plot(real(y2fft));
    set(gca, 'FontSize', 16); title('Real part of FFT');
subplot(222), plot(imag(y1fft)); hold on; plot(imag(y2fft));
    set(gca, 'FontSize', 16); title('Imaginary part of FFT');
subplot(223), plot(abs(y1fft)); hold on; plot(abs(y2fft));
    set(gca, 'FontSize', 16); title('Amplitude');
subplot(224), plot(rad2deg(angle(y1fft))); hold on; plot(rad2deg(angle(y2fft)));
    set(gca, 'FontSize', 16);title('Phase');

% PAUSE   
pause

%%  DEMO 2: demonstrate fftshift with a small vector

%  FFT shift on a vector with an EVEN number of points
test=10*rand(1,6);
testfft=abs(fft(test));         % spit out absolute values of fft
testfftshift=fftshift(testfft);    % unfold the fft
fprintf('FFTSHIFT WITH EVEN NUMBER OF POINTS\n\n');
fprintf('fft amplitudes: \t %.1f %.1f %.1f %.1f %.1f %.1f\n\n', testfft);
fprintf('fftshift contains: \t %.1f %.1f %.1f %.1f %.1f %.1f\n\n', testfftshift);

%  FFT shift on a vector with an ODD number of points
test=10*rand(1,5);
testfft=abs(fft(test));         % spit out absolute values of fft
testfftshift=fftshift(testfft);    % unfold the fft
fprintf('FFTSHIFT WITH ODD NUMBER OF POINTS\n\n');
fprintf('fft amplitudes: \t %.1f %.1f %.1f %.1f %.1f\n\n', testfft);
fprintf('fftshift contains: \t %.1f %.1f %.1f %.1f %.1f\n\n', testfftshift);

% PAUSE
pause

%% DEMO 3:  the back to the cosine function:  unfold the fft and assign wavenumbers

amp2side=abs(fftshift(y1fft))/nx;              % 2-sided FFT amplitudes - have to divide by nx
phs2side=rad2deg(angle(fftshift(y1fft)));      % 2-sided FFT phase in degrees not radians

amp2sidey2=abs(fftshift(y2fft))/nx;              % 2-sided FFT amplitudes - have to divide by nx
phs2sidey2=rad2deg(angle(fftshift(y2fft)));      % 2-sided FFT phase in degrees not radians

figure(3)
subplot(221), plot(amp2side); 
    set(gca, 'FontSize', 16); title('Amplitude');
subplot(222), plot(rad2deg(angle(fftshift(y1fft)))); 
    set(gca, 'FontSize', 16);title('Phase');

% now assign frequencies or wavenumbers
dk=1./(nx*dx);      % freq or wavenumber increments are 1/"period" and "period"=nx*dx;
knyq=1./(2*dx);     % Nyquist wavenumber (or freq, f if time series) = 0.5*sampling freq = 0.5*(1/dx) = 1/(2*dx); 

% THIS IS HOW TO ASSIGN WAVENUMBERS FOR 2-SIDED SPECTRUM - I.E. USE FFT, THEN FFTSHIFT
% AND PLOT BOTH THE POSITIVE AND NEGATIVE WAVENUMBERS)
kk=-knyq:dk:knyq-dk;

subplot(223), plot(kk,amp2side,'k-*'); 
    set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');
subplot(224), plot(kk,rad2deg(angle(fftshift(y1fft)))); 
    set(gca, 'FontSize', 16);title('Phase'); xlabel('Wavenumber (1/m)');

% PAUSE    
pause
    
% NOW make a figure focussing just on the low wavenumbers
figure(4),
subplot(211), plot(kk,amp2side,'k-*'); hold on; plot(kk,amp2sidey2,'r-*');
    set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');
    xlim([-0.04 0.04]); grid on
subplot(212), plot(kk,phs2side,'b-*'); hold on; plot(kk,phs2sidey2,'r-*');
    set(gca, 'FontSize', 16);title('Phase'); xlabel('Wavenumber (1/m)');
    xlim([-0.04 0.04]); grid on

pause
%% DEMO 4:  NOW JUST PLOT THE 1-SIDED (POSITIVE WAVENUMBERS) FFT

% 1-sided FFTs for an EVEN number of points
% take first nx/2+1 points to include k=0 and k=nyquist

amp1side=2*abs(y1fft(1:nx/2+1))/nx;         % twice the amplitude of 2-sided spectrum since combining +ve and -ve freqs  
phs1side=rad2deg(angle(y1fft(1:nx/2+1)));   % phase
kk1=0:dk:knyq;                              % 1-sided wavenumber vector

figure(5),
subplot(211), plot(kk1,amp1side,'k-*'); 
    set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');
    grid on;
subplot(212), plot(kk1,phs1side,'b-*'); 
    set(gca, 'FontSize', 16);title('Phase'); xlabel('Wavenumber (1/m)');
    grid on;

%  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.001);                     % looking for the wavenumber with biggest amplitude
fprintf('Peak amplitude is = %.2f at wavenumber=%.4f \n', amp1side(ik),kk1(ik));

% PAUSE
pause

%% DEMO 5:  FINALLY WHAT IF MY TIME SERIES HAS AN ODD NUMBER OF POINTS - im still not sure this is ok

% clean up
clear all;

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);
amp2side=abs(fftshift(y1fft))/nx;              % 2-sided FFT amplitudes (have to divide by nx)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THIS IS HOW TO ASSIGN WAVENUMBERS 
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-1)*dx);      % wavenumber increment    
kk=-knyq:dk:knyq;       % wavenumbers for 2-sided spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(11),
    subplot(211), plot(kk,amp2side,'k-*'); 
        set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');
        xlim([-0.04 0.04]); grid on
    subplot(212), plot(kk,amp2side,'k-*'); 
        set(gca, 'FontSize', 16); title('Amplitude'); xlabel('Wavenumber (1/m)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  THIS IS HOW TO ASSIGN WAVENUMBERS & select correct # points for 1-sided spectrum
amp1side=2*abs(y1fft(1:(nx+1)/2))/nx;  % 1-sided spectrum: TAKE FIRST (NX+1)/2 POINTS INSTEAD OF (NX/2)+1 POINTS
knyq=1./(2*dx);         % Nyquist wavenumber 
dk=1./(nx*dx);          % wavenumber increment
kk1=0:dk:knyq;          % wavenumbers 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(22),
    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));

% THE END


        