# Symbolic Discrete Fourier Transform (DFT), using FFT?

Is there an implementation of a symbolic FFT in Mathematica?

I've looked around, but I can't find any; Fourier seems to only give numerical results.

• Look at FourierTransform and related functions – Rom38 May 31 '17 at 10:04

The slow ($O(n^2)$) method is to just premultiply your vector with FourierMatrix[]. Using the OP's examples:

FourierMatrix[3, FourierParameters -> {1, -1}].{a, b, c} // ExpToTrig // Simplify
{a + b + c, a - 1/2 I ((-I + Sqrt[3]) b - (I + Sqrt[3]) c),
a + 1/2 I ((I + Sqrt[3]) b - (-I + Sqrt[3]) c)}

FourierMatrix[4, FourierParameters -> {1, -1}].{a, b, c, d} // Simplify
{a + b + c + d, a - I b - c + I d, a - b + c - d, a + I (b + I c - d)}


If you're dealing with small sizes, the difference between $O(n\log n)$ and $O(n^2)$ won't matter too much, so this approach's simplicity might make it preferable in some circumstances.

Note that I had set explicitly set the FourierParameters option to correspond with the expected outputs in the OP's answer. You can of course change it if you have a different preferred normalization.

If the intent is to symbolically analyze the expressions generated in the FFT, then why not just code up a straight implementation at the outset?

symbolicFourier[l_?VectorQ, opts : OptionsPattern[{FourierParameters -> {0, 1}}]] :=
Module[{a, b, n, ω}, {a, b} = OptionValue[FourierParameters];
n = Length[l]; ω = Exp[2 π I b/n];
If[OddQ[n], FourierMatrix[n, opts].l,
Flatten[{{1, 1}, {1, -1}}.{symbolicFourier[l[[1 ;; -2 ;; 2]], opts],
symbolicFourier[l[[2 ;; -1 ;; 2]], opts]
ω^Range[0, Quotient[n, 2] - 1]}]/2^((1 - a)/2)]]


For instance:

symbolicFourier[Array[C, 6], FourierParameters -> {1, -1}]
{C[1] + C[2] + C[3] + C[4] + C[5] + C[6],
C[1] + E^(-((2 I π)/3)) C[3] + E^((2 I π)/3) C[5] +
E^(-((I π)/3)) (C[2] + E^(-((2 I π)/3)) C[4] + E^((2 I π)/3) C[6]),
C[1] + E^((2 I π)/3) C[3] + E^(-((2 I π)/3)) C[5] +
E^(-((2 I π)/3)) (C[2] + E^((2 I π)/3) C[4] + E^(-((2 I π)/3)) C[6]),
C[1] - C[2] + C[3] - C[4] + C[5] - C[6],
C[1] + E^(-((2 I π)/3)) C[3] + E^((2 I π)/3) C[5] -
E^(-((I π)/3)) (C[2] + E^(-((2 I π)/3)) C[4] + E^((2 I π)/3) C[6]),
C[1] + E^((2 I π)/3) C[3] + E^(-((2 I π)/3)) C[5] -
E^(-((2 I π)/3)) (C[2] + E^((2 I π)/3) C[4] + E^(-((2 I π)/3)) C[6])}

• +1 thanks, but it doesn't answer the question. The desire for FFT wasn't related to efficiency... it was so I could get a symbolic expression for the FFT so I could analyze the expression's structure (common subexpressions, etc.)... – user541686 Mar 24 '18 at 3:23
• I guess the answer to your second question is "because this was the only way I knew how to implement an arbitrary-length FFT" (and I'm surprised it's as simple as your code). It's not exactly a trivial problem... – user541686 Mar 24 '18 at 8:36
• I recall you asking about arbitrary-length FFTs in math.SE some time ago, where I did direct you to chirp-z. IIRC, the reason I recommended chirp-z was that I thought you only had a method for power-of-two FFTs. Thankfully there is no such restriction here, but the FFT implementation I gave is still suboptimal, since it does not factor odd factors. If you're interested, I can do the version that factors out the odd part as well, but it might take me longer. – J. M. will be back soon Mar 24 '18 at 8:39
• Oh, I didn't realize it was you who answered that. Yeah, that was quite a long time before this, and indeed for a different purpose, but after that, that was the only way I knew how to do arbitrary-length FFTs, so... yeah. I don't actually need it anymore (the answer I posted here already did what I needed to do.... ironically it got downvoted just an hour-ish ago) so if you post an update it'll have to be for other future readers. :-) – user541686 Mar 24 '18 at 8:50

I have decided to write another answer, since the following method is a vastly different strategy from the one in my previous answer.

The following FFT variant is due to de Boor (of splines fame). Its implementation is pretty compact, at the expense of needing an extra list of the same size as the input. (Contrast this with Cooley-Tukey's classical in-place FFT algorithm.)

deBoorFourier[l_?VectorQ, opts : OptionsPattern[{FourierParameters -> {0, 1}}]] :=
Module[{n = Length[l], a, after, b, before, l1, l2, now, r, ω},
{a, b} = OptionValue[FourierParameters];
l1 = l2 = l; after = 1; before = n;
Do[{now, r} = fac;
Do[before = Quotient[before, now];
Do[ω = Exp[2 π I b (after m + k - 1)/(now after)];
Do[l2[[k + after (m + now j)]] =
Fold[(ω # + #2) &, 0,
Table[l1[[k + after (j + before i)]], {i, now - 1, 0, -1}]],
{j, 0, before - 1}],
{m, 0, now - 1}, {k, after}];
{l1, l2} = {l2, l1}; after *= now,
{r}], {fac, FactorInteger[n]}];
l1/n^((1 - a)/2)]


(Note that I tried to preserve de Boor's variable names in his original FORTRAN code.)

Using the same example in my previous answer:

deBoorFourier[Array[C, 6], FourierParameters -> {1, -1}]
{C[1] + C[2] + C[3] + C[4] + C[5] + C[6],
C[1] - C[4] + E^(-((I π)/3)) (C[2] - C[5] + E^(-((I π)/3)) (C[3] - C[6])),
C[1] + C[4] + E^(-((2 I π)/3)) (C[2] + C[5] + E^(-((2 I π)/3)) (C[3] + C[6])),
C[1] - C[2] + C[3] - C[4] + C[5] - C[6],
C[1] + C[4] + E^((2 I π)/3) (C[2] + C[5] + E^((2 I π)/3) (C[3] + C[6])),
C[1] - C[4] + E^((I π)/3) (C[2] - C[5] + E^((I π)/3) (C[3] - C[6]))}


It is instructive to compare the result of this and my previous routine symbolicFourier[], if at least for its ability to fully split the sequence length into prime factors.

There is a similar FFT method due to Glassman, whose implementation and comparison with de Boor's approach I'll leave to someone else. See also this paper.

• +1 very cool!!! – user541686 Mar 29 '18 at 19:39

...so to be honest, I posted the question as an excuse so I could post this code here. =P

I'm not actually sure if an FFT would be faster symbolically than a direct matrix multiplication, but here's my implementation for 1D FFTs.

Note that the Chirp-Z transform I use for non-power-of-2 FFTs is indirectly based on GNU Octave's, so there may be licensing issues there. I'm not sure if an algorithm like this can be licensed or written any other way, but in any case, if you're redistributing the code, evaluate it and proceed accordingly.

AnyFFT[v_, phasemul_: +1] :=
With[{n = Length[v]},
If[n <= 1, v,
With[{m = n, a = 1},
With[{n2 = If[m + n > 0, 2^BitLength[m + n - 2], 0]},
If[n == 2^(BitLength[n] - 1),
(* FFT code below for power-of-2 size *)
With[
{e = AnyFFT[v[[1 ;; ;; 2]], phasemul],
o = MapIndexed[
{vk, k} \[Function] Exp[-phasemul 2 (Last[k] - 1) I π/n] vk,
AnyFFT[v[[2 ;; ;; 2]], phasemul]]},
(* Chirp Z-Transform below. NOTE: based on GNU Octave's czt.m *)
With[{w = Exp[-phasemul 2 I π/m]},
With[{chirp = w^(Range[1 - n, Max[m, n] - 1]^2/2)},
With[
{xp =
Join[a^-Range[0, n - 1] chirp[[n ;; n + n - 1]] v,
ConstantArray[0, n2 - n]],
icp = Join[1/chirp[[1 ;; m + n - 1]],
ConstantArray[0, n2 - (m + n - 1)]]},
With[{r =
AnyFFT[AnyFFT[xp, phasemul] AnyFFT[icp,
phasemul], -phasemul]/n2},
r[[n ;; m + n - 1]]*chirp[[n ;; m + n - 1]]]]]]]]]]];
AnyIFFT[v_, phasemul_: -1] := AnyFourier[v, -phasemul]/Length[v];


Sample Usage:

AnyFFT[{a, b, c}] // Simplify
AnyFFT[{a, b, c, d}] // Simplify


Sample Output:

\begin{align*} \left\{a+b+c,a-\sqrt[3]{-1} b+(-1)^{2/3} c,a+(-1)^{2/3} b-\sqrt[3]{-1} c\right\} \\ \{a+b+c+d,a-i b-c+i d,a-b+c-d,a+i (b-d)-c\} \end{align*}