7
$\begingroup$

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?

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – J. M. is away Apr 7 '16 at 0:20
8
$\begingroup$

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.

$\endgroup$
  • 2
    $\begingroup$ Just in case this is too clear and easy to read, you can make it obfuscated via dft=With[{n = Length@#, a = #}, (a.Exp[2 π I Range[0, n - 1] #/n]/n & /@ Range[0, n - 1])] & $\endgroup$ – Jason B. Apr 7 '16 at 12:43
  • $\begingroup$ You can compare the result of dftMat with the built-in FourierMatrix[] as another check. $\endgroup$ – J. M. is away Apr 8 '16 at 2:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.