# 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. 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 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 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... 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. 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 Jan 19 '18 at 16:29