# Programming a recursive FFT with Mathematica

I'm trying to program a FFT using a recursive function in Mathematica, though my program is not getting me anywhere at the moment. Could you say what's wrong with it and what I can do to make it work?

My program is the following: (We assume that a is a list of size $2^s$ for some integer $s$)

 FFT[a_] :=

(n := Length[a];
If[n == 1, Return[a],
f := Exp[2 I Pi/n];
o := 1;
b := Drop[a, {1, Length[a], 2}];
c := Drop[a, {2, Length[a], 2}];
d := FFT[b];
e := FFT[c];
For[k = 0, k < n/2, k++,
y[k] := d[[{k}]] + o e[[{k}]];
y[k + n/2] := d[[{k}]] - o e[[{k}]];
o = o f]
Return[Table[y[i], {i, 1, n}]]])

• While I can see there are some issues with your program, I was wondering if you had thought about using the built-in function Fourier instead? – Michael E2 Jul 31 '14 at 17:57
• I have, thanks, but it's more of a programming exercise than anything else. – ThomasZ Jul 31 '14 at 17:59
• Did you write this code yourself or did you translate it from some other language? Is this a homework exercise? If you wrote the code yourself, I'd study the basic of Mathematica a bit before embarking on such a project. There are just too many un-Mathematica-like and inefficient constructs here to go through explaining what's wrong with each of them. For a start, learn about the difference between := and =, don't use For (look up Do and Table), don't use Return when it's not needed, and localize variables with Module. Also check your ;s. – Szabolcs Jul 31 '14 at 19:36
• Thanks for the advice. I wrote it from what I understood mathematically, that's probably why it seems so artificial. That's not homework, I'm just trying to get better at programming with Mathematica, but I guess I'll spend more time on the basics. Cheers, T. – ThomasZ Jul 31 '14 at 21:21
• Also, check your brackets, they are not balanced – Rojo Jul 31 '14 at 21:22

Since you are using this as an excuse to learn to code in Mathematica, I'll try to help with that in mind. Sometimes I find it nice to design "from the top down", document from the top down, split in many small functions with no state, and go testing them "from the bottom up". As you gain confidence, you will perhaps use coarser functions, and not test every little thing.

I'll improvise some "testing framework" for this answer. It's just what came to mind, nothing particularly recommended.

Say we have a function f, and we know that f[2] must return 3, and f[4] must return 9. Then we write test[f] ^= Inactivate@{f[2] -> 3, f[4] -> 9};, and when we want to run those tests, we do runTest[f]. We just need to define something like

runTest[sym_] :=
Inactive[TestReport][
test[sym] /. (x_ -> y_) :> Inactive[VerificationTest][x, y]] //
Activate


Ok, let's start.

So, what do we want to do? An fft. Let's say the fft that's equivalent to Mathematica's Fourier with FourierParameters->{1, -1}.

mmafft=Fourier[#, FourierParameters->{1, -1}]&;


Tests serve as documentation too. Let's write some fft pairs

test[fft] ^= Inactivate@{
fft[{1}] -> {1},
fft[{1, 0, 0, 0}] -> {1, 1, 1, 1},
fft[{0, 0, 0, 1, 0, 0}] -> {1, -1, 1, -1, 1, -1}
};


The fft of a list with a single element is an identity function. Let's write that

fft[{v_}]:={v};


Done. Now, if the list is longer, a generalization of your approach, divide and conquer, would be to split it, apply the fft to each part, and merge the results. Let's write that

fft[l_List]:=fft/@split[l]//merge


Done. Now, we need to tell Mathematica what split and merge are. Let's start with split.

A possible way of splitting, that is consistent with your proposal, is to split "downsampling" to the least integer factor of the list length. So, if the list has 2^s elements, this will always be 2. Let's write some examples

test[split] ^=
Inactivate@{
split[{1}] -> {{1}},
split[{1, 2}] -> {{1}, {2}},
split[Range@6] -> {{1, 3, 5}, {2, 4, 6}}
};


Good. Let's code this

split[l_List]:=split[l, nextFactor[l]];


This definition outsourced the task of knowing in how many pieces to split it. Now, the actual splitting

split[l_List, numParts_Integer]:=l~Partition~numParts//Transpose


This is just a way out of many.

Our split definition relies on a function nextFactor that takes a list and returns the number of pieces to split it into. For example

test[nextFactor] ^= Inactivate@{
nextFactor[Range@10] -> 2,
nextFactor[Range[3 7 2]] -> 2,
nextFactor[Range[3 7]] -> 3};


Let's code this

nextFactor[l_List] := FactorInteger[Length@l][[1, 1]]


We've hit bottom. nextFactor is totally defined. Let's test it. runTest[nextFactor] gives success. Good.

Now, split is defined. Let's test it. runTest[split] gives success. Good.

Time for merge. This function takes a list of ffts from time shifted and downsampled signals and merges them adding them up.

merge[ffts_List] :=
upsampledFft[Length@ffts] /@ ffts // timeShiftedFfts // Total


The fft of an upsampled signal is simply repeating it that many times

test[upsampledFft] ^= Inactivate@{
upsampledFft[2]@{1, 2, 3} -> {1, 2, 3, 1, 2, 3},
upsampledFft[3]@{1, 5} -> {1, 5, 1, 5, 1, 5}
};

upsampledFft[k_Integer][l_List] :=


runTest[upsampledFft] suggests its all good

timeShiftedFfts time shifts each fft a different amount

timeShiftedFfts[ffts_List] := MapIndexed[timeShiftedFft, ffts];


Shifting an fft is multiplying by an exponential that turns as many times in the whole length of the list as the number of shifts.

timeShiftedFft[fft_List, {index_}] := timeShiftedFft[fft, index - 1];
timeShiftedFft[fft_List, n_] := With[{len = Length[fft]},
fft Exp[-I 2 Pi n/len Range[0, len - 1]]
]


I didn't make tests here because I was lazy and they weren't too straightforward, but now we should at least test it with some inputs and see if the output makes sense, starting from timeShfitedFft, up to merge.

merge[{mmafft[{1, 3, 5}], mmafft[{2, 4, 6}]}] // Chop gives the same as mmafft[Range@6], so far so good.

All is done. runTest[fft]. Successsss. This can be made more efficient, but it usually is a good idea to not get crazy about that until you know you need it. Complexity seems fine. Compared to a naive implementation.

naiveFft[l_List] :=
Length@l /. len_ :> l.SparseArray[{i_, j_} :>
Exp[-I 2 Pi (j - 1) (i - 1)/len], len {1, 1}]

r = RandomReal[{-10, 10}, 2^16];
BenchmarkPlot[{naiveFft@# &, fft@# &}, Take[r, #] &,
PowerRange[2^2, 2^16, 2], "IncludeFits" -> True, PlotRange -> Full]


• +1 for showing us the new debugging tools. – Michael E2 Jul 31 '14 at 23:18
• +1. Really nice explanation. Plot added. – RunnyKine Aug 1 '14 at 7:09
• @RunnyKine thanks for the plot :) – Rojo Aug 1 '14 at 7:23
• @Rojo, You're welcome and thanks for the detailed answer. – RunnyKine Aug 1 '14 at 7:31
• @Rojo, might be a dumb question, but what does ^= do? – Sascha Aug 1 '14 at 7:44