7

So to gain a better understanding of how Discrete Fourier Transforms are done, I'm trying to program one in Mathematica so I can easily check my answers with the built-in Mathematica code. Based on Wolfram's page on the DFT, I see that I have to sum over N elements N times to calculate the full transform of a set with N discretizations. So I made a for() loop to calculate the DFT. 

For[k = 0, k < m , k++,
  For[x = 0, x < m , x++,
  Ts[[k + 1]] += s[[x + 1]]*Exp[-2*I*π*k*x/m];
  ]
]

Ts being the transform of s and s being the set that I would like to transform. But every time I check the answer, it's wrong. I thought that I might have a few off-by-one errors but none of the places I would expect there to be an error seem to have one and shifting by one in most places either doesn't make sense or does not product the correct DFT. What am I misunderstanding?

shrx
  • 7,807
  • 2
  • 22
  • 55
Mumbo
  • 71
  • 1
  • 3
    A few points: 1. the default normalization of Fourier[] is not the conventional one; look up FourierParameters. 2. You can use FourierMatrix[] as another possible check of your DFT. 3. Sum[] is built-in, as is Table[]. They should be better than your current loop. 4. Matters Computational is a very good reference for this sort of thing. – J. M.'s missing motivation Apr 07 '16 at 00:20

1 Answers1

8

Here's one way:

x = RandomReal[{-1, 1}, 16];
num = Length[x];
fftx = Table[Total[Table[x[[k + 1]] Exp[2 Pi I k n/num], 
                           {k, 0, num - 1}]], {n, 0, num - 1}]/num;

Now you can check the accuracy by taking the difference between fftx and the built in Fourier command.

Norm[fftx - Fourier[x, FourierParameters -> {-1, 1}], 1]

The difference is typically on the order of 10^-16 or so. Thanks to J.M. for a better way of measuring the difference. As JasonB points out, we can write this more concisely:

m = Length[x];
fftx2 = Total[x[[#]] Exp[2 Pi I Range[0, m - 1] (# - 1)/m] & /@ Range[m]]/m;

Another way to represent the DFT is as a matrix. In this form, the DFT is a matrix multiply and it is easy to see the interpretation as the expansion of the signal x in terms of the (complex-valued) sinusoidal basis functions.

m = Length[x];
dftMat = Table[Exp[2 Pi I k n/m], {n, 0, m - 1}, {k, 0, m - 1}];
dftMat.x

As suggested by J.M., this explicit construction of the DFT matrix can be compared to the built-in function FourierMatrix

dftMat == Sqrt[m] FourierMatrix[m]

which returns True.

bill s
  • 68,936
  • 4
  • 101
  • 191