iPhone ускоряет каркас FFT против Matlab FFT
У меня нет большого математического фона, но часть проекта, над которым я работаю, требует БПФ одного вектора. Функция matlab fft(x) работает точно для того, что мне нужно, но после попытки настроить функции fft в Accelerate Framework я получаю совершенно неточные результаты. Если у кого-то есть больше опыта / опыта в Accelerate Framework fft, я мог бы действительно использовать некоторую помощь, чтобы выяснить, что я делаю неправильно. Я основал свою настройку FFT на примере, который я нашел в Google, но не было никаких учебных пособий или чего-либо, что могло бы привести к другим результатам.
РЕДАКТИРОВАТЬ 1: Изменено вокруг некоторых вещей на основе ответов до сих пор. Кажется, что он делает вычисления, но не выводит их в какой-либо степени близко к Matlab.
Это документация для fft для matlab: 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 Accelerate 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
Цель 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
3 ответа
Одна большая проблема: массивы 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 используют разные определения формулы DFT, особенно в отношении коэффициента масштабирования. Похоже, что БПФ MATLAB включает коррекцию масштабирования 1/N, чего нет в большинстве других БПФ.
Вот полный рабочий пример, выход которого соответствует 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 членами, как это ожидается в случае БПФ с реальным входом.
Также обратите внимание, что нам нужно было масштабировать БПФ vDSP с коэффициентом 2 - это связано с тем фактом, что это БПФ от реального к сложному, которое основано на N/2-точечном БПФ от сложного к сложному, поэтому выходы масштабируются на N/2, тогда как нормальное FFT будет масштабироваться на N.
Я думаю, что это также может быть проблемой упаковки массивов. Я только что посмотрел на их пример кода, и я вижу, что они продолжают вызывать процедуры преобразования, такие как
vDSP_ctoz
Вот пример ссылки на код от Apple: 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}
В MATLAB похоже, что вы делаете fft из 16 реальных значений {1+0i, 2+0i, 3+0i и т. Д.}}, Тогда как в Accelerate вы делаете fft из 8 комплексных значений {1+2i, 3+4i, 5+6i и т.д...}