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&lt;N/2
        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')
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) --&gt; x[1] in sequence
X_e=x(1:2:N);            % even samples begin at x(1)  --&gt; 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&lt;N/2
        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')