iPhone加速框架FFT与Matlab FFT [英] iPhone Accelerate Framework FFT vs Matlab FFT
问题描述
我没有太多的数学背景，但我正在研究的项目的一部分需要单个向量的FFT。 matlab函数fft（x）可以准确地满足我的需要，但在尝试设置Accelerate Framework fft函数后，我得到了完全不准确的结果。如果有人对Accelerate Framework fft有更多的专业知识/经验，我可以真正使用一些帮助来试图弄清楚我做错了什么。我将我的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 *
第1栏至第4栏
1.3600 -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197i
第5至8栏
-0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331i -0.0800 + 0.0159i
第9至12栏
-0.0800 -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i
专栏s 13到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）; //将数据（1- 16）分成赔率和均衡
vDSP_fft_zrip（fftSetup，& fft_data，1，log2n，kFFTDirection_Forward）; // fft forward
vDSP_fft_zrip（fftSetup，& fft_data，1，log2n，kFFTDirection_Inverse）; // fft逆
vDSP_ztoc（& fft_data，2，（COMPLEX *）ffx，1，nOver2）; //将复数组合回实数
目标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
一个大问题：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 *输入;
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（输入\ n）;
for（i = 0; i< n; ++ i）
{
printf（％d：％8g \ n，i，input [i ]）;
}
vDSP_ctozD（（DSPDoubleComplex *）输入，2，& fft_data，1，nOver2）;
printf（FFT Input\\\
）;
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\）;
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\\\
）;
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）; //奈奎斯特
返回0;
}
输出为：
输入
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输入
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输出
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
缩放FFT输出
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
未包装输出
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 =
第1列到第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
第8至14栏：
-8.0000 + 1.5913i -8.0000 + 0.0000i -8.0000 - 1.5913 i -8.0000 - 3.3137i -8.0000 - 5.3454i -8.0000 - 8.0000i -8.0000 - 11.9728i
第15和16栏：
-8.0000 - 19.3137i -8.0000 - 40.2187i
octave-3.4.0：17>
请注意，9到16的输出只是一个复共轭镜像或底8个项，这是实际输入FFT的预期情况。
另请注意，我们需要将vDSP FFT缩放2倍 - 这是因为它是一个真实到复数的FFT，它基于N / 2点复数到复数的FFT，因此输出按N / 2进行缩放，而普通FFT则按N缩放。
I do not have much math background but part of the project I am working on requires the FFT of a single vector. The matlab function fft(x) works accurately for what I need, but after trying to set up the Accelerate Framework fft functions I get completely inaccurate results. If anyone has more expertise/experience with Accelerate Framework fft I could really use some help trying to figure out what I am doing wrong. I based my fft set-up off an example I found on google, but there were no tutorials or anything that produced different results.
EDIT1: Changed around some stuff based on the answers so far. It seems to be doing calculations but it doesnt output them in any way close to that of matlab
This is the documentation for fft for matlab: http://www.mathworks.com/help/techdoc/ref/fft.html
** NOTE: for example purposes, the x array will be {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} in both examples
Matlab Code:
x = fft(x)
Matlab output:
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
Objective C code:
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 output:
ffx now contains:
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
One big problem: C arrays are indexed from 0, unlike MATLAB arrays which are 1-based. So you need to change your loop from
for(int i = 1; i <= 16; i++)
to
for(int i = 0; i < 16; i++)
A second, big problem - you're mixing single precision (float) and double precision (double) routines. Your data is double so you should be using vDSP_ctozD, not vDSP_ctoz, and vDSP_fft_zripD rather than vDSP_fft_zrip, etc.
Another thing to watch out for: different FFT implementations use different definitions of the DFT formula, particularly in regard to scaling factor. It looks like the MATLAB FFT includes a 1/N scaling correction, which most other FFTs do not.
Here is a complete working example whose output matches Octave (MATLAB clone):
#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;
}
Output is:
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
Comparing with Octave we get:
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>
Note that the outputs from 9 to 16 are just a complex conjugate mirror image or the bottom 8 terms, as is the expected case with a real-input FFT.
Note also that we needed to scale the vDSP FFT by a factor of 2 - this is due to the fact that it's a real-to-complex FFT, which is based on an N/2 point complex-to-complex FFT, hence the outputs are scaled by N/2, whereas a normal FFT would be scaled by N.
这篇关于iPhone加速框架FFT与Matlab FFT的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
