Matlab: Radix 2 FFT with 8 Channels and High Number of Inputs

I'm making in Matlab an 8 channel FFT, according to an online tutorial I found at https://www.youtube.com/watch?v=EsJGuI7e_ZQ&t=957s . It works fine for the same number of inputs as shown in the video (N=8), but as I increase N (say for N=128) it no longer works.

I am not sure why it no longer works as I increase N, but I suspect it has something to do with the twiddle factors. I assumed the twiddle factors would have an exponent increasing by 8 after each stage, but maybe they change in a different way as one increases N?

Here is a picture from the YouTube tutorial, at 15:16

Here is a signal flow diagram of my Matlab code

And below is my actual MatLab code. If one uses N=8 (which you can define at top of code) then output is correct. Using N=128 (or other higher values) produces an incorrect result.

I am very confused about what could be causing this problems with the output. Any thoughts?

Much appreciated,

clear all


% Generate input data sequence and plot
N           = 128;
f1          = 10;
num_cycles  = 2;
fs          = f1*N/num_cycles;
x_time      = 0:1/fs:num_cycles/f1-1/fs;
X           = sin(x_time*2*pi*f1);
plot(x_time,X);
title('Input Waveform');


% Split inputs into eight channels
X0_0  =  X(1:8:N);
X0_1  =  X(5:8:N);
X0_2  =  X(3:8:N);
X0_3  =  X(7:8:N);
X0_4  =  X(2:8:N);
X0_5  =  X(6:8:N);
X0_6  =  X(4:8:N);
X0_7  =  X(8:8:N);


% Compute FFT of each channel
X1_0  =  fft(X0_0);
X1_1  =  fft(X0_1);
X1_2  =  fft(X0_2);
X1_3  =  fft(X0_3);
X1_4  =  fft(X0_4);
X1_5  =  fft(X0_5);
X1_6  =  fft(X0_6);
X1_7  =  fft(X0_7);


% Generate Twiddle factors
Wn=exp(-1i*2*pi/N);

% Produce output of first stage of butterfly

for k=0:(N/8)-1
    X2_0(k+1)  =  X1_0(k+1) + (Wn^(k*8)) * X1_1(k+1);
    X2_1(k+1)  =  X1_0(k+1) + (Wn^(k*8+4)) * X1_1(k+1);
    X2_2(k+1)  =  X1_2(k+1) + (Wn^(k*8)) * X1_3(k+1);
    X2_3(k+1)  =  X1_2(k+1) + (Wn^(k*8+4)) * X1_3(k+1);

    X2_4(k+1)  =  X1_4(k+1) + (Wn^(k*8)) * X1_5(k+1);
    X2_5(k+1)  =  X1_4(k+1) + (Wn^(k*8+4)) * X1_5(k+1);
    X2_6(k+1)  =  X1_6(k+1) + (Wn^(k*8)) * X1_7(k+1);
    X2_7(k+1)  =  X1_6(k+1) + (Wn^(k*8+4)) * X1_7(k+1);
end

% Produce output of second stage of butterfly

for k=0:(N/8)-1
    X3_0(k+1)  =  X2_0(k+1) + (Wn^(k*8)) * X2_2(k+1);

    X3_1(k+1)  =  X2_1(k+1) + (Wn^(k*8+2)) * X2_3(k+1);
    X3_2(k+1)  =  X2_0(k+1) + (Wn^(k*8+4)) * X2_2(k+1);
    X3_3(k+1)  =  X2_1(k+1) + (Wn^(k*8+6)) * X2_3(k+1);
    X3_4(k+1)  =  X2_4(k+1) + (Wn^(k*8)) * X2_6(k+1);
    X3_5(k+1)  =  X2_5(k+1) + (Wn^(k*8+2)) * X2_7(k+1);
    X3_6(k+1)  =  X2_4(k+1) + (Wn^(k*8+4)) * X2_6(k+1);
    X3_7(k+1)  =  X2_5(k+1) + (Wn^(k*8+6)) * X2_7(k+1);
end

% Produce output of third stage of butterfly


for k=0:(N/8)-1
    X4_0(k+1)  =  X3_0(k+1) + (Wn^(k*8)) * X3_4(k+1);
    X4_1(k+1)  =  X3_1(k+1) + (Wn^(k*8+1)) * X3_5(k+1);
    X4_2(k+1)  =  X3_2(k+1) + (Wn^(k*8+2)) * X3_6(k+1);
    X4_3(k+1)  =  X3_3(k+1) + (Wn^(k*8+3)) * X3_7(k+1);
    X4_4(k+1)  =  X3_0(k+1) + (Wn^(k*8+4)) * X3_4(k+1);
    X4_5(k+1)  =  X3_1(k+1) + (Wn^(k*8+5)) * X3_5(k+1);
    X4_6(k+1)  =  X3_2(k+1) + (Wn^(k*8+6)) * X3_6(k+1);
    X4_7(k+1)  =  X3_3(k+1) + (Wn^(k*8+7)) * X3_7(k+1);

end

% Merge X4 into final output

c = 1;
for k = 1:N/8
    X5(c:c+7) = [X4_0(k);X4_1(k);X4_2(k);X4_3(k);X4_4(k);X4_5(k);X4_6(k);X4_7(k)];
    c = c+8;
end


% Plot expected FFT and Butterfly FFT
figure
matlab_fft=fft(X);
plot(abs(matlab_fft));
title('Matlab FFT');
figure
plot(abs(X5));
title('Butterfly FFT');

Answer

Your problem is probably with twiddle factor adjustment. Below is a working code for 8 channel implementation of FFT (where the last 3 stages are explicitly computed, but not in a butterfly structure). Its derived based on the recursive structure of FFT. Where the splitting algorithm of N-point DFT into 2 N/2 point DFTs is applied for three consequtive stages.

Mathematical definition of the solution is as follows: Assume that $N$-point DFT $X[k]$ is split into two $N/2$ point DFTS according to the even and odd partitioning as usual: $$ X[k] = X_e[k] + W_N^k X_o[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where the twiddle factor is $W_N = e^{-j\frac{2\pi}{N}}$ and note that the partial DFTs $X_e[k]$ and $X_o[k]$ are of length $N/2$and periodic. Then we can further decompose those $X_e[k]$ and $X_o[k]$ recursively to reach

$$\begin{align} X_e[k] &= X_{ee}[k] + W_{N/2}^k X_{eo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\ X_o[k] &= X_{oe}[k] + W_{N/2}^k X_{oo}[k] ~~~,~~~\text{ for } k=0,1,...,N/2-1 \\ \end{align} $$

And further split those $N/4$ point DFTs into $N/8$ points DFTS which will be the channel signal DFTS.

$$\begin{align} X_{ee}[k] &= X_{eee}[k] + W_{N/4}^k X_{eeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\ X_{eo}[k] &= X_{eoe}[k] + W_{N/4}^k X_{eoo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\ X_{oe}[k] &= X_{oee}[k] + W_{N/4}^k X_{oeo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\ X_{oo}[k] &= X_{ooe}[k] + W_{N/4}^k X_{ooo}[k] ~~~,~~~\text{ for } k=0,1,...,N/4-1 \\ \end{align} $$

Below is an implementation of the above splitting algorithm. Note that your definition of the signals X0_0,X_0_1...,X0_7 is actually correct according to this splitting mode. But the butterfly computations seems not right.

clc;clear all; clear all


% S0 - Generate N point input data sequence x[n] of a sine wave.

% --------------------------------------------------------------
N           = 128;          % Length of signal
M           = 8;            % Number of channels    
x           = sin(2*pi*0.123*[0:N-1]);

% S1 - obtain the even-odd partition signals :
% --------------------------------------------
xe = x(1:2:N);              % EVEN part of N-point x[n]
xo = x(2:2:N);              % ODD  part of N-poibt x[n]


xee  = xe(1:2:N/2);         % EVEN part of N/2-point xe[n]
xeo  = xe(2:2:N/2);         % ODD  part of N/2-point xe[n]
xoe  = xo(1:2:N/2);         % EVEN part of N/2-point xo[n]
xoo  = xo(2:2:N/2);         % ODD  part of N/2-point xo[n]


xeee = xee(1:2:N/4);         % EVEN part of N/4-point xee[n]
xeeo = xee(2:2:N/4);         % ODD  part of N/4-point xee[n]
xeoe = xeo(1:2:N/4);         % EVEN part of N/4-point xeo[n]
xeoo = xeo(2:2:N/4);         % ODD  part of N/4-point xeo[n]

xoee = xoe(1:2:N/4);         % EVEN part of N/4-point xoe[n]
xoeo = xoe(2:2:N/4);         % ODD  part of N/4-point xoe[n]
xooe = xoo(1:2:N/4);         % EVEN part of N/4-point xoo[n]
xooo = xoo(2:2:N/4);         % ODD  part of N/4-point xoo[n]

%S1 - Compute 8 of those N/8 point FFTs for the eee-ooo signals :
% ---------------------------------------------------------------
Xeee = fft(xeee,N/8);
Xeeo = fft(xeeo,N/8);
Xeoe = fft(xeoe,N/8);

Xeoo = fft(xeoo,N/8);
Xoee = fft(xoee,N/8);
Xoeo = fft(xoeo,N/8);
Xooe = fft(xooe,N/8);
Xooo = fft(xooo,N/8);

% S1 - obtain N/4 point DFT's from 8 , N/8 point FFTs of eee-ooo signals :
% ------------------------------------------------------------------------
W4 = exp(-1j*2*pi*(4/N));      % twiddle factor for the 3rd stage (Wn/4)
Xee = [Xeee, Xeee]  + (W4.^[0:N/4-1]).*[Xeeo, Xeeo];

Xeo = [Xeoe, Xeoe]  + (W4.^[0:N/4-1]).*[Xeoo, Xeoo];
Xoe = [Xoee, Xoee]  + (W4.^[0:N/4-1]).*[Xoeo, Xoeo];
Xoo = [Xooe, Xooe]  + (W4.^[0:N/4-1]).*[Xooo, Xooo];

% S2 - obtain N/2 point DFT's from N/4 point FFTs of ee-oo signals :
% ------------------------------------------------------------------
W2 = exp(-1j*2*pi*(2/N));      % twiddle factor for the 2rd stage (Wn/2)
Xe = [Xee, Xee] + (W2.^[0:N/2-1]).*[Xeo, Xeo];
Xo = [Xoe, Xoe] + (W2.^[0:N/2-1]).*[Xoo, Xoo];


% S3 - obtain N point final DFT from two N/2 point DFTs of e-o signals :
% ----------------------------------------------------------------------
Wn = exp(-1j*2*pi/N);      
X = [Xe, Xe] + (Wn.^[0:N-1]).*[Xo, Xo]; 


% SX - Display results
% --------------------
figure,stem(abs(X))
figure,stem(abs(fft(x,N)),'g');