Cooley-Tukey algorithm for Fast Fourier Transform (FFT)

 

Here's a step-by-step explanation of the Cooley-Tukey FFT algorithm:

1. Base Case: If the length of the input sequence x is 1, then the FFT of x is simply x itself. This is the stopping condition for the recursion.

2. Recursive Splitting: If the length of x is greater than 1, the algorithm splits the sequence into two parts: one containing the even-indexed elements and the other containing the odd-indexed elements.

3. FFT Computation: Recursively apply the FFT algorithm to both the even and odd subsequences obtained from the splitting step. This involves repeating steps 1-2 for each subsequence until reaching the base case.

4. Twiddle Factor: Calculate the twiddle factors, which are complex numbers that are applied to the odd-indexed elements before combining them with the even-indexed elements. These factors are calculated as e−2πik/N, where k ranges from 0 to N/2−1, and N is the length of the input sequence.

5. Combine Results: Combine the FFTs of the even and odd subsequences using the twiddle factors. Specifically, multiply the FFT of the odd subsequence by the twiddle factors element-wise, and then add the results to the FFT of the even subsequence. Similarly, subtract the same result from the FFT of the even subsequence.

6. Final Result: The combined results form the FFT of the original sequence xxx. This process continues recursively until the entire FFT is computed.

By recursively applying these steps, the Cooley-Tukey FFT algorithm efficiently computes the FFT of a given sequence.

 

MATLAB Function Code for Cooley-Tukey Fast Fourier Transform (FFT) algorithm

% The code is written by SalimWireless.Com % Cooley-Tukey FFT Algorithm function X = fft_ct(x) N = length(x); if N == 1 X = x; else X_even = fft_ct(x(1:2:end)); % Recursively compute FFT for even-indexed elements X_odd = fft_ct(x(2:2:end)); % Recursively compute FFT for odd-indexed elements factor = exp(-2i * pi * (0:N/2-1) / N); % Twiddle factor X = [X_even + factor .* X_odd, X_even - factor .* X_odd]; % Combine results end end  

 

MATLAB Code for FFT of a Sine Function

% The code is written by SalimWireless.Com clc; clear; close all; % Parameters N = 1024; % Number of samples L = 1; % Length of the interval (in seconds) Fs = 1000; % Sampling frequency (in Hz) f_sine = 50; % Frequency of the sine wave (in Hz) % Time vector t = linspace(0, L, N); % Generate sine wave sine = sin(2 * pi * f_sine * t); % Compute Fourier transform X = fft_ct(sine); % Frequency axis frequencies = linspace(-Fs/2, Fs/2 - Fs/N, N); % Plot plot(frequencies, fftshift(abs(X))); % Shift frequencies for plotting xlabel('Frequency (Hz)'); ylabel('Amplitude'); title('Fourier Transform of Sine Function with 50 Hz'); % Cooley-Tukey FFT Algorithm function X = fft_ct(x) N = length(x); if N == 1 X = x; else X_even = fft_ct(x(1:2:end)); X_odd = fft_ct(x(2:2:end)); factor = exp(-2i * pi * (0:N/2-1) / N); X = [X_even + factor .* X_odd, X_even - factor .* X_odd]; end end  

 Output

 

MATLAB Code for FFT of a Cosine Function

% The code is written by SalimWireless.Com clc; clear; close all; % Parameters N = 1024; % Number of samples L = 1; % Length of the interval (in seconds) Fs = 1000; % Sampling frequency (in Hz) f_sine = 50; % Frequency of the sine wave (in Hz) % Time vector t = linspace(0, L, N); % Generate sine wave sine = cos(2 * pi * f_sine * t); % Compute Fourier transform X = fft_ct(sine); % Frequency axis frequencies = linspace(-Fs/2, Fs/2 - Fs/N, N); % Plot plot(frequencies, fftshift(abs(X))); % Shift frequencies for plotting xlabel('Frequency (Hz)'); ylabel('Amplitude'); title('Fourier Transform of Sine Function with 50 Hz'); % Cooley-Tukey FFT Algorithm function X = fft_ct(x) N = length(x); if N == 1 X = x; else X_even = fft_ct(x(1:2:end)); X_odd = fft_ct(x(2:2:end)); factor = exp(-2i * pi * (0:N/2-1) / N); X = [X_even + factor .* X_odd, X_even - factor .* X_odd]; end end  

 Output

 

MATLAB Code for FFT of a Gaussian Function

% The code is written by SalimWireless.Com clc; clear; close all; % Parameters N = 1024; % Number of samples sigma = 1; % Standard deviation of Gaussian x = linspace(-10, 10, N); % Define x-axis % Gaussian function gaussian = exp(-x.^2 / (2 * sigma^2)); % Compute Fourier transform X = fft_ct(gaussian); % Plot plot(x, abs(X)); xlabel('Frequency'); ylabel('Amplitude'); title('Fourier Transform of Gaussian Function'); % Cooley-Tukey FFT Algorithm function X = fft_ct(x) N = length(x); if N == 1 X = x; else X_even = fft_ct(x(1:2:end)); X_odd = fft_ct(x(2:2:end)); factor = exp(-2i * pi * (0:N/2-1) / N); X = [X_even + factor .* X_odd, X_even - factor .* X_odd]; end end  

Output