8

我没有太多的数学背景，但我正在从事的项目的一部分需要单个向量的 FFT。matlab 函数 fft(x) 可以准确地满足我的需要，但是在尝试设置 Accelerate Framework fft 函数后，我得到的结果完全不准确。如果有人对 Accelerate Framework fft 有更多的专业知识/经验，我真的可以使用一些帮助来找出我做错了什么。我的 fft 设置基于我在 google 上找到的一个示例，但没有教程或任何产生不同结果的东西。

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

苹果加速框架：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

一个大问题：C 数组从 0 开始索引，这与基于 1 的 MATLAB 数组不同。所以你需要改变你的循环

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

到

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

第二个大问题 - 您正在混合单精度 ( float) 和双精度 ( double) 例程。你的数据是double你应该使用vDSP_ctozD的，不是vDSP_ctoz，vDSP_fft_zripD而不是vDSP_fft_zrip，等等。

另一件需要注意的事情：不同的 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 的输出只是复共轭镜像或底部 8 项，与实输入 FFT 的预期情况一样。

另请注意，我们需要将 vDSP FFT 缩放 2 倍——这是因为它是基于 N/2 点复数 FFT 的实数到复数 FFT，因此输出按 N/2 缩放，而正常的 FFT 将按 N 缩放。

于 2012-05-30T19:43:31.960 回答
4

我认为这也可能是数组包装问题。我刚刚查看了他们的示例代码，我看到他们一直在调用转换例程，例如

vDSP_ctoz

这是来自苹果的代码示例链接：http: //developer.apple.com/library/ios/#documentation/Performance/Conceptual/vDSP_Programming_Guide/SampleCode/SampleCode.html#//apple_ref/doc/uid/TP40005147-CH205-中国国际金融学会

我认为这还不是完整的答案，但我也同意 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

In MATLAB, it looks like you're doing an fft of 16 real values {1+0i, 2+0i, 3+0i, etc...} whereas in Accelerate you're doing an fft of 8 complex values {1+2i, 3+4i, 5+6i, etc...}

于 2012-05-31T23:56:31.417 回答