0
$\begingroup$

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. Signal Flow

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

enter image description here

And Matlab's own FFT is shown below

enter image description here

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')
$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Marcus Müller Jan 19 '18 at 10:48
0
$\begingroup$

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<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')
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Karl Haebler Jan 19 '18 at 11:03
  • $\begingroup$ 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 $\endgroup$ – Karl Haebler Jan 19 '18 at 11:08
  • $\begingroup$ 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... $\endgroup$ – Fat32 Jan 19 '18 at 12:17
  • $\begingroup$ 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[1] actually an even sample? I'm getting confused haha. If x[1] is an even sample then you're right. I was under the impression that it's an odd sample. $\endgroup$ – Karl Haebler Jan 19 '18 at 12:27
  • $\begingroup$ @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[0] and even $\endgroup$ – Fat32 Jan 19 '18 at 16:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.