8

私は数学のバックグラウンドがあまりありませんが、私が取り組んでいるプロジェクトの一部では、単一ベクトルの FFT が必要です。matlab 関数 fft(x) は必要なものに対して正確に機能しますが、Accelerate Framework fft 関数をセットアップしようとすると、完全に不正確な結果が得られます。誰かが Accelerate Framework fft の専門知識/経験を持っている場合、私が間違っていることを理解しようとする助けを実際に利用できます。Google で見つけた例に基づいて fft のセットアップを行いましたが、チュートリアルや異なる結果をもたらすものはありませんでした。

EDIT1:これまでの回答に基づいて、いくつかのものを変更しました。計算を行っているようですが、matlab に近い方法でそれらを出力しません

これは、matlab の fft のドキュメントです: http://www.mathworks.com/help/techdoc/ref/fft.html

** 注: 例として、x 配列は両方の例で {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} になります。

Matlab コード:

x = fft(x)

Matlab の出力:

x =

   1.0e+02 *

  Columns 1 through 4

   1.3600            -0.0800 + 0.4022i  -0.0800 + 0.1931i  -0.0800 + 0.1197i

  Columns 5 through 8

  -0.0800 + 0.0800i  -0.0800 + 0.0535i  -0.0800 + 0.0331i  -0.0800 + 0.0159i

  Columns 9 through 12

  -0.0800            -0.0800 - 0.0159i  -0.0800 - 0.0331i  -0.0800 - 0.0535i

  Columns 13 through 16

  -0.0800 - 0.0800i  -0.0800 - 0.1197i  -0.0800 - 0.1931i  -0.0800 - 0.4022i

Apple アクセラレート フレームワーク: http://developer.apple.com/library/mac/#documentation/Accelerate/Reference/vDSPRef/Reference/reference.html#//apple_ref/doc/uid/TP40009464

目的の C コード:

int log2n = log2f(16);

FFTSetupD fftSetup = vDSP_create_fftsetupD (log2n, kFFTRadix2);     

DSPDoubleSplitComplex fft_data;
fft_data.realp = (double *)malloc(8 * sizeof(double));
fft_data.imagp = (double *)malloc(8 * sizeof(double));

vDSP_ctoz((COMPLEX *) ffx, 2, &fft_data, 1, nOver2); //split data (1- 16) into odds and evens

vDSP_fft_zrip (fftSetup, &fft_data, 1, log2n, kFFTDirection_Forward); //fft forward

vDSP_fft_zrip (fftSetup, &fft_data, 1, log2n, kFFTDirection_Inverse); //fft inverse

vDSP_ztoc(&fft_data, 2, (COMPLEX *) ffx, 1, nOver2); //combine complex back into real numbers

客観的な C 出力:

ffx には以下が含まれるようになりました。

272.000000
-16.000000
-16.000000
-16.000000
0.000000
0.000000
0.000000
0.000000
0.000000
10.000000
11.000000
12.000000
13.000000
14.000000
15.000000
16.000000
4

3 に答える 3

15

1 つの大きな問題: C 配列は、1 ベースの MATLAB 配列とは異なり、0 からインデックス付けされます。したがって、ループを次から変更する必要があります

for(int i = 1; i <= 16; i++)

に

for(int i = 0; i < 16; i++)

2 番目の大きな問題 -単精度 ( ) ルーチンfloatと倍精度 ( double) ルーチンを混在させています。あなたのデータはdoubleそうではなく、ではなく、などを使用する必要がありvDSP_ctozDます。vDSP_ctozvDSP_fft_zripDvDSP_fft_zrip

注意すべきもう 1 つのこと: 異なる FFT 実装では、特にスケーリング係数に関して、異なる DFT 式の定義が使用されます。MATLAB FFT には 1/N スケーリング補正が含まれているようですが、他のほとんどの FFT には含まれていません。

出力が Octave (MATLAB クローン) と一致する完全な作業例を次に示します。

#include <stdio.h>
#include <stdlib.h>
#include <Accelerate/Accelerate.h>

int main(void)
{
    const int log2n = 4;
    const int n = 1 << log2n;
    const int nOver2 = n / 2;

    FFTSetupD fftSetup = vDSP_create_fftsetupD (log2n, kFFTRadix2);

    double *input;

    DSPDoubleSplitComplex fft_data;

    int i;

    input = malloc(n * sizeof(double));
    fft_data.realp = malloc(nOver2 * sizeof(double));
    fft_data.imagp = malloc(nOver2 * sizeof(double));

    for (i = 0; i < n; ++i)
    {
       input[i] = (double)(i + 1);
    }

    printf("Input\n");

    for (i = 0; i < n; ++i)
    {
        printf("%d: %8g\n", i, input[i]);
    }

    vDSP_ctozD((DSPDoubleComplex *)input, 2, &fft_data, 1, nOver2);

    printf("FFT Input\n");

    for (i = 0; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }

    vDSP_fft_zripD (fftSetup, &fft_data, 1, log2n, kFFTDirection_Forward);

    printf("FFT output\n");

    for (i = 0; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }

    for (i = 0; i < nOver2; ++i)
    {
        fft_data.realp[i] *= 0.5;
        fft_data.imagp[i] *= 0.5;
    }

    printf("Scaled FFT output\n");

    for (i = 0; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }

    printf("Unpacked output\n");

    printf("%d: %8g%8g\n", 0, fft_data.realp[0], 0.0); // DC
    for (i = 1; i < nOver2; ++i)
    {
        printf("%d: %8g%8g\n", i, fft_data.realp[i], fft_data.imagp[i]);
    }
    printf("%d: %8g%8g\n", nOver2, fft_data.imagp[0], 0.0); // Nyquist

    return 0;
}

出力は次のとおりです。

Input
0:        1
1:        2
2:        3
3:        4
4:        5
5:        6
6:        7
7:        8
8:        9
9:       10
10:       11
11:       12
12:       13
13:       14
14:       15
15:       16
FFT Input
0:        1       2
1:        3       4
2:        5       6
3:        7       8
4:        9      10
5:       11      12
6:       13      14
7:       15      16
FFT output
0:      272     -16
1:      -16 80.4374
2:      -16 38.6274
3:      -16 23.9457
4:      -16      16
5:      -16 10.6909
6:      -16 6.62742
7:      -16  3.1826
Scaled FFT output
0:      136      -8
1:       -8 40.2187
2:       -8 19.3137
3:       -8 11.9728
4:       -8       8
5:       -8 5.34543
6:       -8 3.31371
7:       -8  1.5913
Unpacked output
0:      136       0
1:       -8 40.2187
2:       -8 19.3137
3:       -8 11.9728
4:       -8       8
5:       -8 5.34543
6:       -8 3.31371
7:       -8  1.5913
8:       -8       0

Octave と比較すると、次のようになります。

octave-3.4.0:15> x = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 ]
x =

    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15   16

octave-3.4.0:16> fft(x)
ans =

 Columns 1 through 7:

   136.0000 +   0.0000i    -8.0000 +  40.2187i    -8.0000 +  19.3137i    -8.0000 +  11.9728i    -8.0000 +   8.0000i    -8.0000 +   5.3454i    -8.0000 +   3.3137i

 Columns 8 through 14:

    -8.0000 +   1.5913i    -8.0000 +   0.0000i    -8.0000 -   1.5913i    -8.0000 -   3.3137i    -8.0000 -   5.3454i    -8.0000 -   8.0000i    -8.0000 -  11.9728i

 Columns 15 and 16:

    -8.0000 -  19.3137i    -8.0000 -  40.2187i

octave-3.4.0:17>

9 から 16 までの出力は、実数入力 FFT で予想される場合と同様に、複素共役ミラー イメージまたは下位 8 項であることに注意してください。

また、vDSP FFT を 2 倍にスケーリングする必要があることにも注意してください。これは、N/2 ポイントの複素数から複素数への FFT に基づく実数から複素数への FFT であるためです。出力は N/2 でスケーリングされますが、通常の FFT は N でスケーリングされます。

于 2012-05-30T19:43:31.960 に答える
4

配列のパッキングの問題でもあると思います。彼らのサンプル コードを見たところ、次のような変換ルーチンを呼び出し続けていることがわかりました。

vDSP_ctoz

Appleからのコード サンプル リンクは次のとおりです。 CIAEJIGF

まだ完全な答えではないと思いますが、Paul Rにも同意します.

ちなみに、不思議なことに、Wolfram Alpha に行くと、FFT{1,2,3,4,5,6,7,8,9,10,11,12,13,14, 15,16}

于 2012-05-30T20:09:46.850 に答える
0

MATLAB では、16 の実数値 {1+0i、2+0i、3+0i など...} の fft を実行しているように見えますが、Accelerate では 8 つの複素数値 {1+ の fft を実行しています。 2i、3+4i、5+6i など...}

于 2012-05-31T23:56:31.417 に答える