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<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')
• start as elementary as possible. Test with a zero-vector. If the output is zero, basic sanity. Then start with a vector all zero but for the first entry. is it a constant in the DFT? If not, investigate where you're not adding up correctly. If that works, put that single non-zero entry into the next bin, and so on. Does linearity apply to your implementation? Where does it not? Such questions help you debug stuff. – Marcus Müller Jan 19 '18 at 10:48

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 in sequence
X_e=x(1:2:N);            % even samples begin at x(1)  --> x 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<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')
• Thanks a lot! You're totally right that I should parameterize Wn so that it includes N instead of 32. However, I'm shocked since you changed the sign of the imaginary number to be positive and it worked, whereas from what I read everyone the imaginary number should be negative. How did you know to make the exponent positive? And are you sure that is correct? Thank you very much – Karl Haebler Jan 19 '18 at 11:03
• Hmmmm I see that as one increases the number of cycles (like you did) the effect of the sign change diminishes. But for a low number of cycles (like 2) it is still noticeable. Any thoughts? Thanks – Karl Haebler Jan 19 '18 at 11:08
• First of all of course it is true that $W_N = e^{-j 2 \pi / N}$ hence it was a typo. Then the error encountered after correcting the typo ise due to the mistake of selecting odd and even samples wrongly! Look: x(1:2:N) selects the even samples and x(2:2:N) selects the odd ones... – Fat32 Jan 19 '18 at 12:17
• Hmmmm, but consider the following code A=[1:10]; A_o=A(1:2:10); A_e=A(2:2:10); This properly puts the odds and evens in the right vector, correct? And I'm pretty sure I'm doing the same thing in my main code. Or am I missing something? Or is x actually an even sample? I'm getting confused haha. If x is an even sample then you're right. I was under the impression that it's an odd sample. – Karl Haebler Jan 19 '18 at 12:27
• @KarlHaebler DSP algorithms are mathematically expressed using sequence notation according to which a signal x[n] is a sequence for integer n including ...-2,-1,0,1,2,... Even samples of the sequence is determined by even index n such as -2,0,2,4,... And an dd samples are given by n = ...-3-1,1,3,5.. However the sequence x[n] is represented in MATLAB using a matrix whose indexing start from k=1 always. so x(1) is first sample of the sequence x[n] for its valid samples. Ex. if x[n] is defined for n=0,1,...,N-1 then matrix x(k) will jave index k=1,2,...,N. and x(1) will be x and even – Fat32 Jan 19 '18 at 16:29