% try_fft.m - example use of FFT
% HJSIII - 99.07.28

clear

% create synthetic data
% 2048 samples at 1000 Hz
n = 2048;
fs = 1000;                             % units [Hz]

% 30 Hz sine, +/- 5 V
f = 30;                                % units [Hz]
V = 5;                                 % units [volts]

% optional
% n=2048 and fs=1000 creates incomplete final period
%    spectral resolution misses 30 Hz narrow peak, FFT peaks ~ 3.5 V
%    cleaner phase transitions
% n=2048 and fs=1024 creates integer number of periods
%    spectral resolution hits 30 Hz exactly, FFT peak = 5 V
%    noiser phase transitions
%fs = 1024;                              % units [Hz]

% time
dt = 1 / fs;                            % units [sec]
t = [ 0:(n-1) ]' * dt;                  % units [sec]

% synthetic sine wave
DC = 0;
x = DC + V * sin( 2 * pi * f * t );     % units [volts]

% optional
% white noise with LP filter
%x = V * randn(n,1);
%x = ( x([ 2:n 1 ]) + 2*x + x([ n 1:(n-1) ]) ) / 4; 

% optional
% square wave
%x = DC + V * sign(x);                  % units [volts]

% optional
% 100 mV RMS noise
noise = 0.100;                         % units [volts]
%x = x + noise * randn( size(x) );      % units [volts]

% FFT
% MATLAB FFT scaled by 2/n - DC component scaled by 1/n
a = fft(x) * 2 / n;                    % complex number - units [volts]
a(1) = a(1) / 2;                       % units [volts]
amp = abs( a );                        % units [volts]
phase = angle( a ) * 180 / pi;         % units [deg]
df = fs / n;                           % units [Hz]
freq = [ 0:(n-1) ]' * df;              % units [Hz]

% show results
disp( ' ' )
disp( 'FFT' )
disp( ' freq [Hz] amplitude' )
[ y,i1 ] = max(amp);
i1 = i1 - 5;
i2 = i1 + 10;
disp( [ freq(i1:i2) amp(i1:i2) ] )     % units [Hz] [volts]

% plot time domain
figure( 1 )
subplot( 3,1,1 )
  plot( t,x )                          % units [sec] [volts]
  axis( [ 0  t(n)  -2*V  2*V ] ) 
  xlabel( 'Time (sec)' )
  ylabel( 'Voltage' )
  title( [ num2str(f) ' Hz, +/- ' num2str(V) 'V, ' ...
           num2str(n) ' samples at ' num2str(fs) ' Hz' ] )

% plot amplitude
subplot( 3,1,2 )
  plot( freq,amp )                     % units [Hz] [volts]
  axis( [ 0  freq(n)  0  2*V ] )
  xlabel( 'Frequency [Hz]' )
  ylabel( 'Amplitude [V]' )

% plot phase
subplot( 3,1,3 )
  plot( freq,phase )                   % units [Hz] [deg]
  axis( [ 0  freq(n)  -180  180 ] )
  xlabel( 'Frequency [Hz]' )
  ylabel( 'Phase [deg]' )


% power spectrum AND power spectral density
ipwr = 0;
if ipwr==1,

figure( 2 )

% power spectrum = amplitude squared / 2
ps = amp .* amp / 2;                    % units [volts^2]
% this method is the same
% ps = 4 * fft(x) .* conj(fft(x)) / n / n / 2;
subplot( 3,1,1 )
  plot( freq, ps )                      % units [Hz] [volts^2]
  xlabel( 'Frequency (Hz)' )
  ylabel( 'FFT^2 / 2' )

  disp( ' ' )
  disp( 'FFT^2 / 2' )
  disp( ' freq [Hz]        ps' )
  [ y,i1 ] = max(ps);
  i1 = i1 - 5;
  i2 = i1 + 10;
  disp( [ freq(i1:i2) ps(i1:i2) ] )     % units [Hz] [volts^2]

% power spectral density = power spectrum / spectral bandwidth
%   for equally spaced spectral lines from FFT, spectral 
%     bandwidth = df = fs/n = constant
psd_fft = ps / df;                      % units [volts^2.sec]
subplot( 3,1,2 )
  plot( freq, psd_fft )                 % units [Hz] [volts^2.sec]
  xlabel( 'Frequency [Hz]' )
  ylabel( 'FFT^2 / 2 / df' )

  disp( ' ' )
  disp( 'FFT^2 / 2 / df' )
  disp( ' freq [Hz]   psd_fft' )
  [ y,i1 ] = max(psd_fft);
  i1 = i1 - 5;
  i2 = i1 + 10;
  disp( [ freq(i1:i2) psd_fft(i1:i2) ] )    % units [Hz] [volts^2.sec]

% MATLAB PSD function
% [Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP)
nfft = 256;
[ px,fx ] = psd( x, nfft, fs, nfft, 0 );

subplot( 3,1,3 )
  plot( fx, px )
  xlabel( 'Frequency [Hz]' )
  ylabel( 'MATLAB PSD output' )

  disp( ' ' )
  disp( 'MATLAB PSD output' )
  disp( ' freq [Hz]       psd' )
  [ y,i1 ] = max(px);
  i1 = i1 - 5;
  i2 = i1 + 10;
  disp( [ fx(i1:i2) px(i1:i2) ] )

% bottom of iwpr segment
end