Symbolic Discrete Fourier Transform (DFT), using FFT?

Question

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.

Answer 1

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])}

Answer 2

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.

Answer 3

...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]]},
    Join[MapThread[Plus, {e, o}], MapThread[Subtract, {e, o}]]],
       (* 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*}