Matlab produces two unknown spikes in custom FFT

I'm making a relatively small Cooley–Tukey FFT in Matlab and I'm noticing unusual spikes in the result compared with Matlab's own FFT.

The figure below shows the signal flow of my program. It's a standard Cooley-Tukey scheme.

My results when computing the FFT of a 16 Hz sinusoid are shown below

And Matlab's own FFT is shown below

Clearly I'm missing something. What could be the cause of the two large spikes in the middle? I conjecture that it has something to do with how the even and odd parts are combined, for instance if there is a discontinuity there. I'm really not sure though.

Any help is greatly appreciated!

The code I'm using is as follows

clear all

% Generate input data sequence and plot
N=128;
f1=16;

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

% split inputs into even and odd samples and compute fft of each division
X_o=x(1:2:N);
X_e=x(2:2:N);


fft_x_o=fft(X_o);
fft_x_e=fft(X_e);

% Generate base twiddle factor
W32=exp(-1i*2*pi/32);

% Combine fft even and odd with twiddle factors to produce final output
for k=0:N-1
    if k        X(k+1)=fft_x_e(k+1)+(W32^k)*fft_x_o(k+1);

    else
        X(k+1)=fft_x_e(k+1-N/2)+(W32^k)*fft_x_o(k+1-N/2);
    end
end

% plot butterfly fft and matlab fft
FFT_xaxis=0:fs/N:fs-fs/N;
figure
plot(FFT_xaxis,abs(X))
title('Butterfly FFT')

xlabel('Frequency')
ylabel('Magnitude')
matlab_fft=fft(x);
figure
plot(FFT_xaxis,abs(matlab_fft))
title('Matlab FFT')
xlabel('Frequency')
ylabel('Magnitude')

Answer

The following is the corrected code. It seems you have the problem in the twiddle factor and the selection of even and odd samples of x[n]...

N=128;
f1=16;
num_cycles=32;
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);

% split inputs into even and odd samples and compute fft of each division
X_o=x(2:2:N);            % odd samples begin at x(2) --> x[1] in sequence

X_e=x(1:2:N);            % even samples begin at x(1)  --> x[0] in sequence

fft_x_o = fft(X_o, N/2);
fft_x_e = fft(X_e, N/2);

% Generate base twiddle factor
WN=exp(-1i*2*pi/N);

% Combine fft even and odd with twiddle factors to produce final output
for k=0:N-1

    if k        X(k+1) = fft_x_e(k+1)+(WN^k)*fft_x_o(k+1);
    else
        X(k+1) = fft_x_e(k+1-N/2)+(W3N^k)*fft_x_o(k+1-N/2);
    end
end

% plot butterfly fft and matlab fft
FFT_xaxis=0:fs/N:fs-fs/N;
figure

plot(FFT_xaxis,abs(X))
title('Butterfly FFT')
xlabel('Frequency')
ylabel('Magnitude')
matlab_fft=fft(x);
figure
plot(FFT_xaxis,abs(matlab_fft))
title('Matlab FFT')
xlabel('Frequency')
ylabel('Magnitude')