N13page 1 of 12
Fourier Transforms
Discrete Fourier Transform (DFT)
find frequency content inherent within a signal
guess specific frequencies fk and assess how well each matches the given signal
digitize xi for i = 1 to n at fixed sampling frequency fs = 1/h
choose frequency fk
create guess gi for i = 1 to n using frequency fk and compare to xi
correlation coefficient
correlation coefficient for sine at frequency fk
Ck is approximate amplitude at frequency fk inherent within digitized signal xi
correlation coefficient
convolution integral over time sample T
unknownsignal at frequency fo
guess
perfect guess fk = fo
time T for m = 1, 2, 3, … complete cycles at frequency fo
form cycles plus (or minus) some small portion of another cycle
DFT at guess fk requires correlation for both sine and cosine
DC fk = 0 Ak = 2 mean(xi) Bk = 0
1) single DFT at one given fk
2) multiple DFT at several discrete frequencies over band f1 < fk < f2
3) complete DFT at n discrete frequencies for k = 0 to n
spectral resolution
small number of samples n causes coarse resolution in frequency domain
why n discrete frequencies ?
Order of Computational Complexity
FLOP = floating point operation = one multiply and one add
assume sine and cosine functions have been pre-computed for given number of samples n
O[ ] = number of FLOPS for computation
O[ Ak ] = n O[ Bk ] = n O[ single DFT ] = 2n
O[ complete DFT ] = 2n(n+1) ≈ 2n2
larger sample size n provides finer spectral resolution
larger sample size n increases computation time by n2
Fast Fourier Transform (FFT)
Cooley and Tukey (1965)
radix-2algorithm breaks n point DFT into two smaller (n/2) point DFTs
compute DFT of even index points ( x2 x4 x6 … )
compute DFT of odd index points ( x1 x3 x5 … )
combine results using symmetric “butterfly” operations
recursively break each (n/2) DFTs into smaller (n/4) DFTs
must use n = 2P for P = integer
n = …, 1024, 2048, 4096, …
provides significantly faster performance over DFT
O[ complete DFT ] = 2 n2 n = 1024 requires 2,097,152 FLOPS
O[ FFT ] = n ln(n) n = 1024 requires 7,098 FLOPS
What about 3000 data points?
1) discard samples at beginning or end and use n = 2048
2) pad extra values onto the end and use n = 4096
a) pad with zeros
b) pad with last value
3) use non-radix-2 FFT (slower)
Practical considerations of using FFT
1) not all FFT libraries normalize by 2/n
2) most FFT libraries use complex values real = Ak imaginary = Bk
3) DFT and FFT are fully reversible
% try_fft.m - example use of FFT
% HJSIII - 13.11.18
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
% noisier phase transitions
%fs = 1024; % units [Hz]
% time
h = 1 / fs; % units [sec]
t = [ 0:(n-1) ]' * h; % 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]
% optional
% exponential decay
tau = 2.0;
envelope = exp( -3 * t / tau );
%x = x .* envelope;
% 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;
ifipwr==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