FFT recursive code problem

Hi,

I wrote a recursive FFT program. But when I run the program, the result is different from the build-in fft in MATLAB. Can someone help me to find the problems in the code? Thanks!

Jian

function F = fft_rec(f)
n = length(f);  
if (n == 1)
  F = f;
else
  f_even = f(1:2:n);
  f_odd = f(2:2:n);
  X1 = fft_rec(f_even);
  X2 = fft_rec(f_odd).*Wn(n);  
  F1 = X1 + X2;
  F2 = X1 - X2;
  F = [F1 F2];
end

function W = Wn(n)
m = n/2;
W = exp(-2*pi*j.*[0:1:m-1]/n);

function X = my_fft(f)
  n = length(f);
  N = pow2(ceil(log2(n)));
  f = [f, zeros(1, N - n)];
  X = fft_rec(f);
  X = X(1:n);

function test_fft()
  t = 1 : 40;
  y = sin(t);
  F = zeros(1, length(t));
  F = my_fft(y);
  x = zeros(1, length(F));
  for i = 1 : length(F)
    x(i) = i;
  end
  plot(x, abs(F));
  hold;
  F2 = fft(y);
  plot(x, abs(F2),'c');
  plot(x, abs(F2)-abs(F), 'r');

Answer 1

[Walter Roberson]

I'm not sure what X2 gets multiplied by Wn(n) but X1 does not?

If n is always odd for Wn(n) then why are you confusing things by using m = n/2 knowing that will be a number with a fraction, and then trying to use that fraction as the upper boundary of the 0:1:m-1 array?

Answer 2

[Jian]

Hi Walter,

Thank you for your response. For the Wn(n), the n is not odd, it is the length of the sequence. Actually, it should be 2^n. For recursive fft, it divide the sequence into even and odd parts, then calculate each part, and compose the result.

Best,

Jian

[Walter Roberson]

You have if (n==1) and do something there. The implication is that it is possible for n==3 or other odd values. When that odd value (say 7) is passed to Wn(n) then m = 7/2 = 3.5, w = exp(-2*pi*j.*[0:1:3.5]/7) which will be exp(-2i*pi.*[0 3]/7) which will be [exp(0) exp(-6i/7*pi). f_odd will be f([2 4 6]) in this situation but f_even will be f([1 3 5 7]) in this situation, which is going to lead to crashes when you go to add the vector of length 4 to the vector of length 3.

Unless, that is, you posted incorrect code.

[Jian]

Hi Walter,

Good point. I made the sequence whose length is 2^n and solved the problem recursively. The subsequence lengths becomes ..16, 8, 4, 2, 1. For n == 1, I do calculate and return. The Wn only deals with 2, 4, 8, 16, ... . Does it make sense? Actually, I referred some DSP books, like "Digital Signal and Processing Principal, Algorithms, and Applications>> , <<Introduction to Algorithms>>(MIT Press). I checked my code, but can not get the hint. If you have any clue, please let me know.

Thank you.

Answer 3

[David Young]

The return variable for the function Wn is W, but the result is assigned to w. Make them both the same, and it works fine.

Ideally, there should be a check that the length of the input is even (if it is not 1), to avoid a mysterious error message if the original length is not a power of 2.

(Please also format code in questions using the {}Code button so that the lines don't get run together.)

[Jian]

Hi David,

Thank you for your response. In my code, for "Wn" function, it assigns to "W" and returns. I typed wrong "w" in this post. Sorry for that. It looks like it is not the error of assignment.

If you have any clue, please point it out and help me to solve the problem.

Also, thanks for the advice to use {}. I am new to here and didn't know that.

Jian

[David Young]

Hi Jian,

When I changed "W" to "w", your code worked perfectly. That is, it gave the same answer as the built-in fft function, to within a tolerance of about 10^-15.

What goes wrong for you - is there an error message, or do you get a numerical difference?

Thanks for editing the question to make the code clear.

[Jian]

Hi David,

Thank you for your response and testing my code. I really appreciate that. In order to test my code, I used sin wave, for instance, x = [0:100], y = sin(x), and compare the result with my_fft(y) and fft(y), (which fft is build-in function in MATLAB). The numerical is different. I also did some pad to make sure the input sequence always being 2^n. In other words, you don't need to worry about the unbalanced sequences. The result is not satisfied as you said. I edited my original post and pasted the test code. You can check it out if you are free. I guess the problem may reside in the sampling. But I can not and have not testify it in math. BTW, I used the student version MATLAB.

Best

[David Young]

The difference you see is because MY_FFT pads the input to 2^n entries, but FFT does not. If you change t=1:40 in your test code to t=1:128 (or any other power of 2), so no padding is needed, you'll find you get the results you expect. Alternatively, pad the input signal before you pass it to FFT.

[Jian]

Hi David,

I followed your advice and it worked. Thanks! To handle the input sequence whose length is not 2^n, I used 0 to pad the sequence. But I guess MATLAB build-in fft cuts the sequence to be 2^n instead of padding the sequence with 0. That is why it can always get the satisfied result. Does my guess make sense?

Jian

[Walter Roberson]

No, Matlab does not use 2^n sequences: it factors the length of the array and then works on fragments determined by the factorization. I do not know the details.

[Jian]

Hi Walter,

Would you mind explain the "factors" and "fragment" more? For example, in order to handle a sequence with length 40, does MATLAB cut the sequence to be length 32, or does it estimate some data and make the sequence to be length 48?

Thanks!

[Walter Roberson]

Neither. It is not necessary for to use powers of two in order to fft. It might, for example, factor it in to 5 * 8, use the power-of-2 algorithm on the 8 portion, and use a different fft algorithm for the 5 part.

[Jian]

Thank you for the explanation. It is very helpful. As my understanding, the radix-2 fft is mostly used and very efficient. That is why I assumed that MATLAB uses it.

[Walter Roberson]

http://www.dspguru.com/dsp/faqs/fft

1.6 Are FFTs limited to sizes that are powers of 2?

No. The most common and familiar FFTs are "radix 2". However, other radices are sometimes used, which are usually small numbers less than 10. For example, radix-4 is especially attractive because the "twiddle factors" are all 1, -1, j, or -j, which can be applied without any multiplications at all.

Also, "mixed radix" FFTs also can be done on "composite" sizes. In this case, you break a non-prime size down into its prime factors, and do an FFT whose stages use those factors. For example, an FFT of size 1000 might be done in six stages using radices of 2 and 5, since 1000 = 2 * 2 * 2 * 5 * 5 * 5. It might also be done in three stages using radix 10, since 1000 = 10 * 10 * 10.

[Jian]

Thanks!