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

1) discard samples at beginning or end and use n = 2048

2) pad extra values onto the end and use n = 4096

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