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