iPhone Accelerate Framework FFT vs Matlab FFT

Question

У меня нет много математического фона, но часть проекта, над которым я работаю, требует БПФ одного вектора. Функция matlab fft (x) работает точно для того, что мне нужно, но, пытаясь настроить функции FLE ускорения Framework, я получаю совершенно неточные результаты. Если у кого-то есть больше опыта/опыта с Accelerate Framework fft, я мог бы действительно использовать некоторую помощь, пытаясь понять, что я делаю неправильно. Я основывал свою настройку fft на примере, который я нашел в Google, но не было никаких учебников или чего-либо, что давало разные результаты.

EDIT1: Изменено вокруг некоторых вещей, основанных на ответах до сих пор. Это, кажется, делает расчеты, но оно не будет выводить их любым способом, близким к MATLAB

Это документация для FFT для MATLAB: http://www.mathworks.com/help/techdoc/ref/fft.html

** Примечание: для примера, х массив будет быть {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 Ускорить Framework: 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 

Objective выход 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 

Answer 1

Одна большая проблема: C массивы индексируются от 0, в отличие от массивов MATLAB которые основаны на 1. Таким образом, вы должны изменить свой цикл из

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 используют разные определения формулы ДПФ, в частности, в отношении коэффициента масштабирования , Похоже, что MATLAB FFT включает в себя коррекцию масштабирования 1/N, которой нет в большинстве других БПФ.

---

Вот полный рабочий пример, выход которого соответствует октаву (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-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 с реальным входом.

Обратите также внимание на то, что нам нужно было масштабировать VFTSP-FFT в 2 раза - это связано с тем, что это реалистичный комплексный БПФ, основанный на комплексном комплексе N/2 FFT, поэтому выходы масштабируются на N/2, тогда как нормальный БПФ будет масштабироваться N.

Answer 2

Я думаю, что это также может быть проблемой упаковки массива. Я просто смотрел на их пример кода, и я вижу, что они продолжают вызова процедуры преобразования как

vDSP_ctoz 

Вот образец ссылку код из яблока: http://developer.apple.com/library/ios/#documentation/Performance/Conceptual/vDSP_Programming_Guide/SampleCode/SampleCode.html#//apple_ref/doc/uid/TP40005147-CH205-CIAEJIGF

Я не думаю, что это полный ответ все же, но я также согласен с Полом Р.

Кстати, просто любопытно, если вы пойдете в Wolfram Alpha, они дадут совершенно другой ответ для FFT {1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16}

Answer 3

В MATLAB, похоже, что вы делаете n fft из 16 действительных значений {1 + 0i, 2 + 0i, 3 + 0i и т. д.}, тогда как в Accelerate вы выполняете fft из 8 комплексных значений {1 + 2i, 3 + 4i, 5 + 6i, etc ...}