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;
    }

そして実際の計算では、おそらくこれがバグの場所です。

    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引数に含まれます。

ここに画像の説明を入力してください

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フレームワークライブラリで次のコードを見つけました。コードで処理されていない追加の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リンクでは、各値に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 に答える