12

我正在尝试实现离散傅立叶变换，但它不起作用。我可能在某个地方写了一个错误，但我还没有找到它。

基于以下公式：

大地

此函数执行第一个循环，循环 X0 - Xn-1... 在此处输入图像描述

    public Complex[] Transform(Complex[] data, bool reverse)
    {
        var transformed = new Complex[data.Length];
        for(var i = 0; i < data.Length; i++)
        {
            //I create a method to calculate a single value
            transformed[i] = TransformSingle(i, data, reverse);
        }
        return transformed;
    }

而实际计算，这可能就是bug所在。

    private Complex TransformSingle(int k, Complex[] data, bool reverse)
    {
        var sign = reverse ? 1.0: -1.0;
        var transformed = Complex.Zero;
        var argument = sign*2.0*Math.PI*k/data.Length;
        for(var i = 0; i < data.Length; i++)
        {
            transformed += data[i]*Complex.FromPolarCoordinates(1, argument*i);
        }
        return transformed;
    }

接下来解释其余代码：

var sign = reverse ? 1.0: -1.0;反转的 DFT 不会-1在参数中包含，而常规 DFT 确实-1在参数中包含 a。

在此处输入图像描述

var argument = sign*2.0*Math.PI*k/data.Length;是算法的参数。这部分：

在此处输入图像描述

然后是最后一部分

transformed += data[i]*Complex.FromPolarCoordinates(1, argument*i);

我想我仔细复制了算法，所以我看不出我在哪里犯了错误......

附加信息

正如 Adam Gritt 在他的回答中所表明的那样，AForge.net 对该算法有一个很好的实现。我大概可以通过复制他们的代码在 30 秒内解决这个问题。但是，我仍然不知道我在实施中做错了什么。

我真的很好奇我的缺陷在哪里，以及我的解释是错误的。

4

2 回答 2

5

我现在做复杂数学的日子已经过去了，所以我自己可能会遗漏一些东西。但是，在我看来，您正在执行以下操作：

transformed += data[i]*Complex.FromPolarCoordinates(1, argument*i);

什么时候应该更像：

transformed += data[i]*Math.Pow(Math.E, Complex.FromPolarCoordinates(1, argument*i));

除非您将其包含在方法中FromPolarCoordinates()

更新：我在AForge.NET Framework库中找到了以下代码，它显示了正在执行的其他 Cos/Sin 操作，这些操作没有在您的代码中处理。此代码可以在 Sources\Math\FourierTransform.cs: DFT 方法的完整上下文中找到。

for ( int i = 0; i < n; i++ )
{
    dst[i] = Complex.Zero;

    arg = - (int) direction * 2.0 * System.Math.PI * (double) i / (double) n;

    // sum source elements
    for ( int j = 0; j < n; j++ )
    {
        cos = System.Math.Cos( j * arg );
        sin = System.Math.Sin( j * arg );

        dst[i].Re += ( data[j].Re * cos - data[j].Im * sin );
        dst[i].Im += ( data[j].Re * sin + data[j].Im * cos );
    }
}

它使用自定义的 Complex 类（因为它是 4.0 之前的版本）。大多数数学与您实现的相似，但内部迭代正在对实部和虚部进行额外的数学运算。

进一步更新：经过一些实施和测试，我发现上面的代码和问题中提供的代码产生了相同的结果。根据评论，我还发现此代码生成的内容与 WolframAlpha 生成的内容之间有什么区别。结果的不同之处在于，Wolfram 似乎对结果应用了 1/sqrt(N) 的归一化。在所提供的 Wolfram Link 中，如果每个值都乘以 Sqrt(2)，则这些值与上述代码生成的值相同（舍入误差除外）. 我通过将 3、4 和 5 个值传递给 Wolfram 对此进行了测试，发现我的结果分别与 Sqrt(3)、Sqrt(4) 和 Sqrt(5) 不同。基于离散傅里叶变换维基百科提供的信息确实提到了使 DFT 和 IDFT 的变换成为单一的规范化。这可能是您需要考虑修改代码或了解 Wolfram 可能在做什么的途径。

于 2011-04-19T12:10:37.277 回答
1

您的代码实际上几乎是正确的（您在逆变换中缺少 1/N）。问题是，您使用的公式通常用于计算，因为它更轻，但在纯理论环境（以及在 Wolfram 中），您将使用 1/sqrt(N) 的归一化来使变换成为单一的。

即你的公式是：

Xk = 1/sqrt(N) * sum(x[n] * exp(-i*2*pi/N * k*n))

x[n] = 1/sqrt(N) * sum(Xk * exp(i*2*pi/N * k*n))

这只是归一化的约定问题，只有幅度会发生变化，因此您的结果还不错（如果您没有忘记逆变换中的 1/N）。

干杯

于 2011-04-20T21:21:35.350 回答