Fast Fourier transform (FFT) MATLAB tutorial series (Part 2.1) - - PowerPoint PPT Presentation

▶

fast fourier transform fft

Fast Fourier transform (FFT) MATLAB tutorial series (Part 2.1) - - PowerPoint PPT Presentation

Aug 15, 2022 177 likes •709 views

Fast Fourier transform (FFT) MATLAB tutorial series (Part 2.1) Pouyan Ebrahimbabaie Laboratory for Signal and Image Exploitation (INTELSIG) Dept. of Electrical Engineering and Computer Science University of Lige Lige, Belgium Applied

slide-1

SLIDE 1

Fast Fourier transform (FFT)

Pouyan Ebrahimbabaie Laboratory for Signal and Image Exploitation (INTELSIG)

Dept. of Electrical Engineering and Computer Science

University of Liège Liège, Belgium Applied digital signal processing (ELEN0071-1) 1 April 2020

MATLAB tutorial series (Part 2.1)

slide-2

SLIDE 2

Introduction 𝒖 |𝒀 𝒈 | 𝒈 𝒚(𝒖)

slide-3

SLIDE 3

Introduction 𝒖 |𝒀 𝒈 | 𝒈 𝒚(𝒖)

DFT

slide-4

SLIDE 4

Introduction 𝒖 |𝒀 𝒈 | 𝒈 𝒚(𝒖)

DFT

Time: sampled, finite Frequency: sampled, finite

slide-5

SLIDE 5

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Time: sampled, finite Frequency: sampled, finite

slide-6

SLIDE 6

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Each sample in frequency domain corresponds to a sample in time domain !

slide-7

SLIDE 7

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Time: N samples Frequency: N samples

slide-8

SLIDE 8

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Time: N samples Frequency: N samples

𝑼𝒕

slide-9

SLIDE 9

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Time: N samples Frequency: N samples

𝑼𝒕 𝑮𝒕 = 𝟐 𝑼𝒕

slide-10

SLIDE 10

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Time: N samples Frequency: N samples

𝑼𝒕 𝑮𝒕 = 𝟐 𝑼𝒕 𝑼𝒏𝒃𝒚

slide-11

SLIDE 11

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 | 𝒈

DFT

Time: N samples Frequency: N samples

𝑼𝒕 𝑮𝒕 = 𝟐 𝑼𝒕 𝑼𝒏𝒃𝒚 = 𝑶 − 𝟐 × 𝑼𝒕 𝑼𝒏𝒃𝒚

slide-12

SLIDE 12

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 |

DFT

Time: N samples Frequency: N samples

𝑼𝒕 𝑮𝒕 = 𝟐 𝑼𝒕 𝑼𝒏𝒃𝒚 = 𝑶 − 𝟐 × 𝑼𝒕 𝑼𝒏𝒃𝒚 𝑮𝒕/𝟑 −𝑮𝒕/𝟑

From sampling theorem

slide-13

SLIDE 13

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 |

DFT

Time: N samples Frequency: N samples

𝑼𝒕 𝑮𝒕 = 𝟐 𝑼𝒕 𝑼𝒏𝒃𝒚 = 𝑶 − 𝟐 × 𝑼𝒕 𝑼𝒏𝒃𝒚 𝑮𝒕/𝟑 −𝑮𝒕/𝟑 𝑮𝒕

slide-14

SLIDE 14

Introduction 𝒚(𝒖) 𝒖 |𝒀 𝒈 |

DFT

Time: N samples Frequency: N samples

𝑼𝒕 𝑮𝒕 = 𝟐 𝑼𝒕 𝑼𝒏𝒃𝒚 = 𝑶 − 𝟐 × 𝑼𝒕 𝑼𝒏𝒃𝒚 𝑮𝒕/𝟑 −𝑮𝒕/𝟑 𝑮𝒕

𝑮𝒕/(𝑶 − 𝟐)

slide-15

SLIDE 15

FFT is just an algorithm to compute DFT efficiently ! In MATLAB: X=fft(x)

r

X=fft(x,n)

slide-16

SLIDE 16

Example 2.1: pure tone % Sampling frequency Fs=44100; % Sampling period Ts=1/Fs; % Signal length (in second) N_sec=5; % Siganal length (in sample) N=N_secFs; % Maximum time Tmax=(N-1)Ts; % Time vector t=0:Ts:Tmax;

slide-17

SLIDE 17

Example 2.1: pure tone % Sampling frequency Fs=44100; % Sampling period Ts=1/Fs; % Signal length (in second) N_sec=5; % Siganal length (in sample) N=N_secFs; % Maximum time Tmax=(N-1)Ts; % Time vector t=0:Ts:Tmax;

slide-18

SLIDE 18

Example 2.1: pure tone % Sampling frequency Fs=44100; % Sampling period Ts=1/Fs; % Signal length (in second) N_sec=5; % Siganal length (in sample) N=N_secFs; % Maximum time Tmax=(N-1)Ts; % Time vector t=0:Ts:Tmax;

slide-19

SLIDE 19

Example 2.1: pure tone % Sampling frequency Fs=44100; % Sampling period Ts=1/Fs; % Signal length (in second) N_sec=5; % Siganal length (in sample) N=N_secFs; % Maximum time Tmax=(N-1)Ts; % Time vector t=0:Ts:Tmax;

slide-20

SLIDE 20

Example 2.1: pure tone % Signal x=sin(2piF0.*t); % Play the sound sound(x,Fs) % Plot the sound (show from zero to 60 msec) figure(1) plot(t,x,'LineWidth',2.5) xlim([0 0.06])

slide-21

SLIDE 21

Example 2.1: pure tone % Signal x=sin(2piF0.*t); % Play the sound sound(x,Fs) % Plot the sound (show from zero to 60 msec) figure(1) plot(t,x,'LineWidth',2.5) xlim([0 0.06]) % take FFT (without shift) X1=fft(x);

slide-22

SLIDE 22

Example 2.1: pure tone % Signal x=sin(2piF0.*t); % Play the sound sound(x,Fs) % Plot the sound (show from zero to 60 msec) figure(1) plot(t,x,'LineWidth',2.5) xlim([0 0.06]) % take FFT (without shift) X1=fft(x); % plot result figure(2) plot(abs(X1),'LineWidth',2.5);

slide-23

SLIDE 23

Example 2.1: pure tone % Signal x=sin(2piF0.*t); % Play the sound sound(x,Fs) % Plot the sound (show from zero to 60 msec) figure(1) plot(t,x,'LineWidth',2.5) xlim([0 0.06]) % take FFT (without shift) X1=fft(x); % plot result figure(2) plot(abs(X1),'LineWidth',2.5);

You should always perform some post processing operations (shifting, scaling, etc.) to be able to present the results of fft !

slide-24

SLIDE 24

Example 2.1: pure tone % take FFT (with shift) X2=fftshift(fft(x)); % Frequency range F=-Fs/2:Fs/(N-1):Fs/2; figure(3) plot(F,abs(X2)/N,'LineWidth',2.5); xlabel('Frequency (Hz)') title('Double sided magnitude response')

slide-25

SLIDE 25

Example 2.1: pure tone % take FFT (with shift) X2=fftshift(fft(x)); % Frequency range F=-Fs/2:Fs/(N-1):Fs/2; figure(3) plot(F,abs(X2)/N,'LineWidth',2.5); xlabel('Frequency (Hz)') title('Double sided magnitude response')

slide-26

SLIDE 26

Example 2.1: pure tone % take FFT (with shift) X2=fftshift(fft(x)); % Frequency range F=-Fs/2:Fs/(N-1):Fs/2; figure(3) plot(F,abs(X2)/N,'LineWidth',2.5); xlabel('Frequency (Hz)') title('Double sided magnitude response')

slide-27

SLIDE 27

Example 2.1: pure tone % take FFT (with shift) X2=fftshift(fft(x)); % Frequency range F=-Fs/2:Fs/(N-1):Fs/2; figure(3) plot(F,abs(X2)/N,'LineWidth',2.5); xlabel('Frequency (Hz)') title('Double sided magnitude response') fftshift Scaling

slide-28

SLIDE 28

Example 2.1: pure tone % take FFT (with shift) X2=fftshift(fft(x)); % Frequency range F=-Fs/2:Fs/(N-1):Fs/2; figure(3) plot(F,abs(X2)/N,'LineWidth',2.5); xlabel('Frequency (Hz)') title('Double sided magnitude response')

This is the correct method to graph a “double-sided” (negative and positive) frequency spectrum

slide-29

SLIDE 29

Example 2.2: single sided spectrum % Sampling frequency Fs=44100; % Sampling period Ts=1/Fs; % Signal length (in second) N_sec=5; % Signal length (in sample) N=N_secFs; % Maximum time Tmax=(N-1)Ts; % Time vector t=0:Ts:Tmax;

slide-30

SLIDE 30

Example 2.2: single sided spectrum % pure F0, F1 and F3 F0=600; F1=1300; F2=2000; % Signal x=sin(2piF0.t)+… …0.5sin(2piF1.t)+0.2sin(2piF2.*t); % Play the sound sound(x,Fs) % Plot the sound (show from zero to 60 msec) figure(1) plot(t,x,'LineWidth',2.5) xlim([0 0.06])

slide-31

SLIDE 31

Example 2.2: single sided spectrum % Compute fft X=fft(x); % Take abs and scale it X2=abs(X/N); % Pick the first half X1=X2(1:N/2+1); % Multiply by 2 (except the DC part), to compensate % the removed side from the spectrum. X1(2:end-1) = 2X1(2:end-1); % Frequency range F = Fs(0:(N/2))/N;

slide-32

SLIDE 32

Example 2.2: single sided spectrum % Compute fft X=fft(x); % Take abs and scale it X2=abs(X/N); % Pick the first half X1=X2(1:N/2+1); % Multiply by 2 (except the DC part), to compensate % the removed side from the spectrum. X1(2:end-1) = 2X1(2:end-1); % Frequency range F = Fs(0:(N/2))/N;

slide-33

SLIDE 33

Example 2.2: single sided spectrum % Compute fft X=fft(x); % Take abs and scale it X2=abs(X/N); % Pick the first half X1=X2(1:N/2+1); % Multiply by 2 (except the DC part), to compensate % the removed side from the spectrum. X1(2:end-1) = 2X1(2:end-1); % Frequency range F = Fs(0:(N/2))/N;

slide-34

SLIDE 34

Example 2.2: single sided spectrum % Compute fft X=fft(x); % Take abs and scale it X2=abs(X/N); % Pick the first half X1=X2(1:N/2+1); % Multiply by 2 (except the DC part), to compensate % the removed side from the spectrum. X1(2:end-1) = 2X1(2:end-1); % Frequency range F = Fs(0:(N/2))/N;

slide-35

SLIDE 35

Example 2.2: single sided spectrum % Compute fft X=fft(x); % Take abs and scale it X2=abs(X/N); % Pick the first half X1=X2(1:N/2+1); % Multiply by 2 (except the DC part), to compensate % the removed side from the spectrum. X1(2:end-1) = 2X1(2:end-1); % Frequency range F = Fs(0:(N/2))/N;

slide-36

SLIDE 36

Example 2.2: single sided spectrum % Plot single-sided spectrum plot(F,X1,'LineWidth',2.5) title('Single-Sided Amplitude Spectrum') xlabel('f (Hz)')

Most of the time we use “single-sided” amplitude or phase spectrum !

slide-37

SLIDE 37

Example 2.3: siren F0=1300; F1=200; F2=1400; B0=100; B1=100; B2=500; % Signal x=sin(2piF0.t+B0pit.^2)+… sin(2piF1.t+B1pit.^2)... +sin(2piF2.t+B2pi*t.^2);

slide-38

SLIDE 38

Example 2.4: voice % read audio file .wav [x,Fs]=audioread('adult_female_speech.wav'); % play the sound sound(x,Fs) % Sampling period Ts=1/Fs; % Length of signal N=length(x); % Maximum time Tmax=(N-1)*Ts; % Time vector t=0:Ts:Tmax; …

slide-39

SLIDE 39

Usable voice frequency band in telephony: ~ 300 Hz to 3400 Hz

slide-40

SLIDE 40

How fast is it?

𝒀(𝒇𝒌𝝏) = ෍

𝒐=𝟏 𝑶−𝟐

𝒚[𝒐]𝒇−𝒌𝝏𝒐 DTFT:

slide-41

SLIDE 41

How fast is it?

𝒀(𝒇𝒌𝝏) = ෍

𝒐=𝟏 𝑶−𝟐

𝒚[𝒐]𝒇−𝒌𝝏𝒐 DTFT: Sample 𝝏 at 𝝏𝒍 = Τ (𝟑𝝆 𝑶)𝒍, 𝒍 = 𝟏, 𝟐, 𝟑, ⋯ 𝑶

slide-42

SLIDE 42

How fast is it?

𝒀(𝒇𝒌𝝏) = ෍

𝒐=𝟏 𝑶−𝟐

𝒚[𝒐]𝒇−𝒌𝝏𝒐 DTFT: Sample 𝝏 at 𝝏𝒍 = Τ (𝟑𝝆 𝑶)𝒍, 𝒍 = 𝟏, 𝟐, 𝟑, ⋯ 𝑶 𝒀 𝒍 = 𝒀(𝒇𝒌𝟑𝝆

𝑶 𝒍) = ෍ 𝒐=𝟏 𝑶−𝟐

𝒚[𝒐]𝒇−𝒌𝟑𝝆

𝑶 𝒍𝒐

DFT:

slide-43

SLIDE 43

How fast is it?

𝒀(𝒇𝒌𝝏) = ෍

𝒐=𝟏 𝑶−𝟐

𝒚[𝒐]𝒇−𝒌𝝏𝒐 DTFT: Sample 𝝏 at 𝝏𝒍 = Τ (𝟑𝝆 𝑶)𝒍, 𝒍 = 𝟏, 𝟐, 𝟑, ⋯ 𝑶 𝒀 𝒍 = 𝒀(𝒇𝒌𝟑𝝆

𝑶 𝒍) = ෍ 𝒐=𝟏 𝑶−𝟐

𝒚[𝒐]𝒇−𝒌𝟑𝝆

𝑶 𝒍𝒐

DFT: For each 𝒍: 𝑶 complex multiplications, 𝑶 − 𝟐 complex adds

slide-44

SLIDE 44

How fast is it?

𝑷(𝑶𝟑) computations for direct DFT

slide-45

SLIDE 45

How fast is it?

𝑷(𝑶𝟑) computations for direct DFT vs. 𝑷(𝑶 𝐦𝐩𝐡𝟑 𝑶) computations for fft !

slide-46

SLIDE 46

How fast is it?

𝑷(𝑶𝟑) computations for direct DFT vs. 𝑷(𝑶 𝐦𝐩𝐡𝟑 𝑶) computations for fft ! 𝑶 𝟐𝟏𝟏𝟏 𝟐𝟏𝟕 𝟐𝟏𝟘 𝑶𝟑 𝟐𝟏𝟕 𝟐𝟏𝟐𝟑 𝟐𝟏𝟐𝟗 𝑶 𝐦𝐩𝐡𝟑 𝑶 𝟐𝟏𝟓 𝟑𝟏 × 𝟐𝟏𝟕 𝟒𝟏 × 𝟐𝟏𝟘

slide-47

SLIDE 47

How fast is it?

𝑷(𝑶𝟑) computations for direct DFT vs. 𝑷(𝑶 𝐦𝐩𝐡𝟑 𝑶) computations for fft ! 𝑶 𝟐𝟏𝟏𝟏 𝟐𝟏𝟕 𝟐𝟏𝟘 𝑶𝟑 𝟐𝟏𝟕 𝟐𝟏𝟐𝟑 𝟐𝟏𝟐𝟗 𝑶 𝐦𝐩𝐡𝟑 𝑶 𝟐𝟏𝟓 𝟑𝟏 × 𝟐𝟏𝟕 𝟒𝟏 × 𝟐𝟏𝟘 𝟐𝟏𝟐𝟗 𝒐𝒕 ~ 𝟒𝟐. 𝟑 𝒛𝒇𝒃𝒔𝒕 vs. 𝟒𝟏 × 𝟐𝟏𝟘 𝒐𝒕 ~ 𝟒𝟏 𝒕𝒇𝒅𝒑𝒐𝒆𝒕

slide-48

SLIDE 48

First time presented by …

Cooley and Tukey (1965)

slide-49

SLIDE 49

First time presented by …

Gauss (1805)

slide-50

SLIDE 50

Real-Time application …

slide-51

SLIDE 51

Useful links

https://www.youtube.com/watch?v=iTMn0Kt18tg
https://allsignalprocessing.com/fast-fourier-transform-

fft-algorithm/

https://allsignalprocessing.com/discrete-fourier-

transform-sampling-the-dtft/

https://nl.mathworks.com/help/matlab/ref/fft.html