Understanding the radix-2 FFT recursive algorithm

Question

I'm currently studying DFT and FFT and we were given this simple question: se the recursive FFT to compute the DFT of this polynomial of third degree: -1 + 4x + 3x^2.

So, I'm considering this algorithm:

How does this recursion work? Is the for loop positioned at the end of all the recursion calls? Or every time y goes back to y^0 and y^1? Can someone guide me? Of course, not all the steps, just some examples? Thank you so much in advance!

Answer 1

The recursive implementation of the radix-2 Decimation In Frequency algorithm can be understood using the following two figures. The first one refers to pushing the stack phase, while the second one illustrates the popping the stack phase.

In particular, the two figures illustrate the following Matlab implementation:

% --- Radix-2 Decimation In Frequency - Iterative approach

function y = radix2_DIF_Recursive(x)

global sumCount mulCount

N = length(x);
phasor = exp(-2 * pi * 1i / N) .^ (0 : N / 2 - 1);

if N == 1
    y = x;
else
    y_top       = radix2_DIF_Recursive(x(1: 2 : (N - 1)));
    y_bottom    = radix2_DIF_Recursive(x(2 : 2 : N));
    z           = phasor .* y_bottom;
    y           = [y_top + z, y_top - z];
    sumCount    = sumCount + 6 * (N / 2);
    mulCount    = mulCount + 4 * (N / 2);
end

Use the following main script to test it:

% --- Radix-2 Decimation In Frequency - Iterative approach

clear all
close all
clc

global sumCount mulCount

N = 32;

x = randn(1, N);

sumCount = 0;
mulCount = 0;

tic
xhat = radix2_DIF_Recursive(x);
timeCooleyTukey = toc;

tic
xhatcheck = fft(x);
timeFFTW = toc;

rms = 100 * sqrt(sum(sum(abs(xhat - xhatcheck).^2)) / sum(sum(abs(xhat).^2)));

fprintf('Time Cooley-Tukey = %f; \t Time FFTW = %f\n\n', timeCooleyTukey, timeFFTW);
fprintf('Theoretical multiplications count \t = %i; \t Actual multiplications count \t = %i\n', ...
         2 * N * log2(N), mulCount);
fprintf('Theoretical additions count \t\t = %i; \t Actual additions count \t\t = %i\n\n', ...
         3 * N * log2(N), sumCount);
fprintf('Root mean square with FFTW implementation = %.10e\n', rms);

The above recursive implementation is the counterpart of the iterative implementation of the radix-2 Decimation In Frequency at Implementation of the Discrete Fourier Transform - FFT. Please, note that the present recursive one does not require any reverse bit ordering, neither of the input nor of the output sequences.