I'm writing a very simple in-place DFT. I am using the formula shown here: http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Definition along with Euler's formula to avoid having to use a complex number class just for this. So far I have this:
private void fft(double[] data)
{
double[] real = new double[256];
double[] imag = new double[256];
double pi_div_128 = -1 * Math.PI / 128;
for (int k = 0; k < 256; k++)
{
for (int n = 0; n < 256; n++)
{
real[k] += data[k] * Math.Cos(pi_div_128 * k * n);
imag[k] += data[k] * Math.Sin(pi_div_128 * k * n);
}
data[k] = Math.Sqrt(real[k] * real[k] + imag[k] * imag[k]);
}
}
But the Math.Cos and Math.Sin terms eventually go both positive and negative, so as I'm adding those terms multiplied with data[k], they cancel out and I just get some obscenely small value. I see how it is happening, but I can't make sense of how my code is perhaps mis-representing the mathematics. Any help is appreciated. FYI, I do have to write my own, I realize I can get off-the shelf FFT's.
I believe You have an error in this section
for (int n = 0; n < 256; n++)
{
real[k] += data[k] * Math.Cos(pi_div_128 * k * n);
imag[k] += data[k] * Math.Sin(pi_div_128 * k * n);
}
You should replace data[k] with data[n]
EDIT:
You are also destroying Your data in the line:
data[k] = Math.Sqrt(real[k] * real[k] + imag[k] * imag[k]);
You have to store the moduli of your complex numbers somewhere else or later. If the modulus is all you want, and You want to store it in data[] You have to code another loop after computing the transform., The whole data[] is necessary for computing every real[k] and imag [k].
private double[] dft(double[] data)
{
int n = data.Length;
int m = n;// I use m = n / 2d;
double[] real = new double[n];
double[] imag = new double[n];
double[] result = new double[m];
double pi_div = 2.0 * Math.PI / n;
for (int w = 0; w < m; w++)
{
double a = w * pi_div;
for (int t = 0; t < n; t++)
{
real[w] += data[t] * Math.Cos(a * t);
imag[w] += data[t] * Math.Sin(a * t);
}
result[w] = Math.Sqrt(real[w] * real[w] + imag[w] * imag[w]) / n;
}
return result;
}
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With