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');

• exactly the same what I commented to your last question: "Incorrect Result" isn't really sufficient description of the problem. Build simple test cases! Don't use us as your debugger! – Marcus Müller Jan 22 '18 at 12:34
• Thanks for the feedback Marcus. To take your advice, I've used two simple test cases to debug. One is an 8 point sin wave. I cannot use the vector with all zeroes except one element like you suggest since the input has to be periodic. So my first test case is a success (for N=8, everything works). But increasing N causes problems, which makes me think my issue is with the twiddle factor. My second test case is in fact the script related to my previous post (which is the same thing I'm doing here except with two channels). That previous problem was so that I understood how to better make this – Karl Haebler Jan 22 '18 at 12:55
• Basically my main question is how do the twiddle factors change after each iteration of 8 inputs. I'm assuming that they are exponentiated by a factor of 8, but I can't find online articles or other sources confirming this. Or more generally what modifications to the FFT are required for a higher N, since I am reading a lot of online articles but they all deal with lower N, and not for iterated inputs with N=8. Does that make sense? – Karl Haebler Jan 22 '18 at 12:56
• Ok basically what I mean is forget about the structure of the Matlab code. My problem is that I am unsure about how the twiddle factors should change for each iteration of 8 inputs... – Karl Haebler Jan 22 '18 at 16:12
• Have you tried different twiddle factor exponent variations? – Fat32 Jan 22 '18 at 17:29

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');

• thanks so much for your help! That's really really appreciated, especially showing the math and your code. After long consideration I've ditched the idea of always keeping data in 8 separate channels. Like you show, I'm now combining the data in half fewer channels after each stage (so starting with 8 channels, then 4, then 2, then 1). For anyone else interested I'm posting the signal flow diagram of this solution, so you can compare it with my original signal flow diagram. I do have one remaining question which I'll ask in the next comment !Signal Flow – Karl Haebler Jan 23 '18 at 14:46
• You're solution makes sense and thanks a bunch for that. One thing I'm curious about though is where you concatenate Xe[k] and Xo[k] with themselves before combining to calculate the output X[k]. I get that Xe[k] and Xo[k] are half as long (length N/2) as X[k] (length N), but how would I know that they are supposed to be concatenated? Do you know what mathematical premise this is described by? Thanks! – Karl Haebler Jan 23 '18 at 15:03
• As I said, this is not a butterfly solution, rather a recursive one, (eventhough the two are inherently related so far). Yes you're right, since the even/odd components of any stage are half long the stage, they should be repeated by two periods to properly sum up the stage. And one way of doing this is to concatenate them... – Fat32 Jan 23 '18 at 15:14