Calculate Discrete Fourier Transform by using c#

There is a definition of Discrete Fourier Transform:

I want to calculate amplitude and phase of a transformed sequence.

I tried :

 public ArrayList phase = new ArrayList();
public ArrayList amp = new ArrayList();
public float phi = 0;

void Start ()
{
DFT();
foreach (var ampList in amp)
{
Debug.Log(ampList);
}
}

void DFT()
{
//Given sequence of Numbers

float[]  X = {1,2,3};
float[]  Y = {2,3,5};
float  N = X.Length;
float re = 0;
float im = 0;

for (int k = 0; k < N; k++)
{
//The discrete Fourier transform, lists of Re & Im instead of a+ib
for (int n = 0; n < N; n++)
{
float phi =  (Mathf.PI * k *n)/N;
re += X[n] * Mathf.Cos(phi) - Y[n] * Mathf.Sin(phi);
im -= Y[n] * Mathf.Cos(phi) - X[n] * Mathf.Sin(phi);
}

re = re/N;
im = im/N;
var  ampp = Mathf.Sqrt(re*re + im*im);
var  phasse  = Mathf.Atan2(re,im);

phase.Add(phasse);
amp.Add(ampp);
}
}


but answer is Not correct, I use Wolfram Mathematica to check the answer :

• This looks like a general programming problem that doesn't require game-specific expertise to answer, since discrete Fourier transforms are used in many non-game fields. You may find you get answers faster on our general programming sister site StackOverflow due to their higher traffic. Or you may find more Fourier experts on the Digital Signal Processing StackExchange. Just be sure to delete the copy here if you choose to ask elsewhere — cross-posting is discouraged on the SE network of sites. – DMGregory Dec 10 '19 at 17:11

1 Answer

There are a number of minor issues that are getting in your way. The following corrections are assuming you are using the definition you provided by Wikipedia:

1. The re and im variables should be inside the for loop for variable k, so that they get reset for every element of the output sequence
2. phi is missing a factor of 2
3. The signs of some of the sin and cos terms are reversed
4. There is no need to divide re and im by N

The final code should look something more like this:

float[] X = {1, 2, 3};
float[] Y = {2, 3, 5};
float N = X.Length;

for (int k = 0; k < N; k++)
{
float re = 0;
float im = 0;

//The discrete Fourier transform, lists of Re & Im instead of a+ib
for (int n = 0; n < N; n++)
{
float phi = (2 * Mathf.PI * k * n) / N;
re += X[n] * Mathf.Cos(phi) + Y[n] * Mathf.Sin(phi);
im += Y[n] * Mathf.Cos(phi) - X[n] * Mathf.Sin(phi);
}

var ampp = Mathf.Sqrt(re * re + im * im);
var phasse = Mathf.Atan2(re, im);

phase.Add(phasse);
amp.Add(ampp);
}


Also, Wolfram Mathematica's Fourier function uses a different convention by default than the one you are using. See the section on FourierParameters in the documentation of the Fourier function. It looks like the definition found in Wikipedia is using the "signal processing" convention, and so you would need to try this in Wolfram Mathematica instead.

Fourier[{1 + 2 * I, 2 + 3 * I, 3 + 5 * I}, FourierParameters -> {1, -1}]


Unfortunately I do not currently have access to Wolfram Mathematica for testing, and Wolfram Alpha does not seem to support this syntax according to this stack exchange, so I will leave that part up to you :-)

As a final note, while it is definitely good to implement the Discrete Fourier Transform yourself as a learning exercise, I would recommend searching for a library which implements this for use in an actual application. Especially if you are going to be processing a lot of data, then a "Fast Fourier Transform" implementation in a good library is likely to be faster, better tested, and more accurate than anything you could make yourself in a reasonable amount of time.

Hope this helps, and best of luck!