Message Boards Message Boards

2 Replies
2 Total Likes
View groups...
Share this post:

Issues while attempting to implement my own version of FFT

Posted 1 year ago

Hey everyone. As the title suggests, I am attempting to implement a manual version of the Cooley-Tukey Fast Fourier Algorithm in Mathematica but I cannot get it to work. Below is the code that I have written:

FFT[list_] := If[Length[list] <= 2, Fourier[list],
  N2 := Length[list];
  Wn := E^(2 \[Pi] I/N2);
  W := 1;
  Aeven := list[[2 ;; ;; 2]];
  Aodd := list[[1 ;; ;; 2]];
  Yeven := FFT[Aeven];
  Yodd := FFT[Aodd];
  For[j = 0, j < N/2, j++,
   Y[j] := Yeven[j] + W*Yodd[j];
   Y[j + N/2] := Yeven[j] - W*Yodd[j];
   W := W*Wn];]]

I am unsure how to get the last For-loop to save the values in the new list Y, and how to then return it so that it can be accessed later. For reference, this is some pseudocode of what I am trying to do: enter image description here

Any help and/or comments would be greatly appreciated!

POSTED BY: Koen Denekamp
2 Replies
Posted 1 year ago

The main problem with your code is that it keeps overwriting the variables during the recursion. You need to localize those variables (I use Module below). Another minor issue is that you can just use Set (=) rather than SetDelayed (:=) in most of your cases. Another problem is that you used N as a variable (I assume it should be N2, since you initialized N2 earlier). The reason why you're not getting any output is that For is a special form that always returns Null. You'll need to explicity return a value. Also, you can simplify the length condition with overloaded DownValues rather than using an If.

I don't have test cases, so there may be mistakes in the details, but here is something that is probably more along the lines of what you want:

newFFT[list : {_}] := list;
newFFT[list_List /; IntegerQ[Log2[Length@list]]] :=
  {Wn = E^(2 \[Pi] I/Length[list]),
   Yeven = newFFT[list[[1 ;; ;; 2]]],
   Yodd = newFFT[list[[2 ;; ;; 2]]],
   Yresult = ConstantArray[0, Length@list]},
   j = 1; W = 1; N2 = Length[list],
   j <= N2/2,
   Yresult[[j]] = Yeven[[j]] + W*Yodd[[j]];
   Yresult[[j + N2/2]] = Yeven[[j]] - W*Yodd[[j]];
   W *= Wn];

One confusing thing here is that I kept the variable names Yeven and Yodd, but Yeven actually is the odd-indexed elements and Yodd is the even ones. In the pseudo code, indexing starts at 0, so Yeven includes the first element. Mathematica indexing starts at 1. I could have switched the names, but I wanted to keep things aligned with the pseudocode. You should probably come up with your own names to clarify the algorithm.

If this doesn't produce the right output (I notice that it doesn't match the built in Fourier function, but I know there are various definitions, so maybe that's okay), then provide some test cases to debug with.

By the way, I'm assuming that you were wanting to implement the algorithm faithfully, but just in case you weren't aware, there are built-in functions Fourier functions. There might also be more idiomatic ways to implement this in Mathematica.

POSTED BY: Eric Rimbey
Posted 1 year ago

Thank you so much! I figured out that the reason why it gives different answers, is that the algorithm does not normalize, while the default Fourier[] command does do some normalization, so adding the "FourierParameters -> {1,1}" option, which removes the normalization, does give the same values (Up to some numerical differences).

POSTED BY: Koen Denekamp
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract