希尔伯特变换（分析信号）使用苹果的加速框架？
我在使用Apple的Accelerate Framework与C++相当的Hilbert变换的问题。 我已经能够得到vDSP的FFTalgorithm，并且在Paul R的post的帮助下，得到了与Matlab相同的结果。
我已经阅读了这两个： 乔丹这个stackoverflow问题，并阅读了Matlabalgorithm（在“algorithm”子标题下） 。 总结algorithm分三个阶段：
- 进行input的FFT。
- 直stream和奈奎斯特之间的零reflection频率和双倍频率。
- 对修改的正向FFT输出进行逆FFT。
下面是每个阶段的Matlab和C ++的输出。 这些示例使用以下数组：
- Matlab：
m = [1 2 3 4 5 6 7 8]; - C ++：
float m[] = {1,2,3,4,5,6,7,8};
Matlab示例
阶段1：
36.0000 + 0.0000i -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 + 0.0000i -4.0000 - 1.6569i -4.0000 - 4.0000i -4.0000 - 9.6569i
阶段2：
36.0000 + 0.0000i -8.0000 + 19.3137i -8.0000 + 8.0000i -8.0000 + 3.3137i -4.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
阶段3：
1.0000 + 3.8284i 2.0000 - 1.0000i 3.0000 - 1.0000i 4.0000 - 1.8284i 5.0000 - 1.8284i 6.0000 - 1.0000i 7.0000 - 1.0000i 8.0000 + 3.8284i
C ++示例（使用Apple的Accelerate框架）
阶段1：
Real: 36.0000, Imag: 0.0000 Real: -4.0000, Imag: 9.6569 Real: -4.0000, Imag: 4.0000 Real: -4.0000, Imag: 1.6569 Real: -4.0000, Imag: 0.0000
阶段2：
Real: 36.0000, Imag: 0.0000 Real: -8.0000, Imag: 19.3137 Real: -8.0000, Imag: 8.0000 Real: -8.0000, Imag: 3.3137 Real: -4.0000, Imag: 0.0000
阶段3：
Real: -2.0000, Imag: -1.0000 Real: 2.0000, Imag: 3.0000 Real: 6.0000, Imag: 7.0000 Real: 10.0000, Imag: 11.0000
很明显，最终的结果是不一样的。 我期望能够计算Matlab的“阶段3”结果（或者至less是虚构的部分），但我不确定如何去做，我尝试了所有我能想到的没有成功的事情。 我感觉这是因为我没有将Apple Accelerate版本中的reflection频率归零，因为它们没有提供（由于是从DC和Nyquist之间的频率计算出来的），所以我认为FFT只是采用共轭加倍的频率作为reflection值，这是错误的。 如果有人能帮我解决这个问题，我将不胜感激。
码
void hilbert(std::vector<float> &data, std::vector<float> &hilbertResult){ // Set variable. dataSize_2 = data.size() * 0.5; // Clear and resize vectors. workspace.clear(); hilbertResult.clear(); workspace.resize(data.size()/2+1, 0.0f); hilbertResult.resize(data.size(), 0.0f); // Copy data into the hilbertResult vector. std::copy(data.begin(), data.end(), hilbertResult.begin()); // Perform forward FFT. fft(hilbertResult, hilbertResult.size(), FFT_FORWARD); // Fill workspace with 1s and 2s (to double frequencies between DC and Nyquist). workspace.at(0) = workspace.at(dataSize_2) = 1.0f; for (i = 1; i < dataSize_2; i++) workspace.at(i) = 2.0f; // Multiply forward FFT output by workspace vector. for (i = 0; i < workspace.size(); i++) { hilbertResult.at(i*2) *= workspace.at(i); hilbertResult.at(i*2+1) *= workspace.at(i); } // Perform inverse FFT. fft(hilbertResult, hilbertResult.size(), FFT_INVERSE); for (i = 0; i < hilbertResult.size()/2; i++) printf("Real: %.4f, Imag: %.4f\n", hilbertResult.at(i*2), hilbertResult.at(i*2+1)); }
上面的代码是我用来获得“Stage 3”C ++（带有Apple的Accelerate Framework）结果的。 前向fft的fft(..)方法执行vDSP fft程序，然后按0.5的比例进行解压缩（按照Paul R的post）。 当执行反转fft时，数据被打包，使用vDSP缩放2.0，最后按1 /（2 * N）缩放。
所以主要问题（就我所知，因为你的代码示例不包括对vDSP的实际调用）是你正在尝试使用从第三步到现实复杂的FFT程序，这基本上是复数到复数的逆FFT（如你的Matlab结果具有非零虚数部分的证据）。
下面是一个使用vDSP的简单的C实现，它匹配你期望的Matlab输出（我使用了更先进的vDSP_DFT例程，这个例程通常应该比较老的fft例程更受欢迎，否则这与你所做的非常相似，对于实数到复数的正向变换，但是从复数到复数的逆变换）：
#include <Accelerate/Accelerate.h> #include <stdio.h> int main(int argc, char *argv[]) { const vDSP_Length n = 8; vDSP_DFT_SetupD forward = vDSP_DFT_zrop_CreateSetupD(NULL, n, vDSP_DFT_FORWARD); vDSP_DFT_SetupD inverse = vDSP_DFT_zop_CreateSetupD(forward, n, vDSP_DFT_INVERSE); // Look like a typo? The real-to-complex DFT takes its input separated into // the even- and odd-indexed elements. Since the real signal is [ 1, 2, 3, ... ], // signal[0] is 1, signal[2] is 3, and so on for the even indices. double even[n/2] = { 1, 3, 5, 7 }; double odd[n/2] = { 2, 4, 6, 8 }; double real[n] = { 0 }; double imag[n] = { 0 }; vDSP_DFT_ExecuteD(forward, even, odd, real, imag); // At this point, we have the forward real-to-complex DFT, which agrees with // MATLAB up to a factor of two. Since we want to double all but DC and NY // as part of the Hilbert transform anyway, I'm not going to bother to // unscale the rest of the frequencies -- they're already the values that // we really want. So we just need to move NY into the "right place", // and scale DC and NY by 0.5. The reflection frequencies are already // zeroed out because the real-to-complex DFT only writes to the first n/2 // elements of real and imag. real[0] *= 0.5; real[n/2] = 0.5*imag[0]; imag[0] = 0.0; printf("Stage 2:\n"); for (int i=0; i<n; ++i) printf("%f%+fi\n", real[i], imag[i]); double hilbert[2*n]; double *hilbertreal = &hilbert[0]; double *hilbertimag = &hilbert[n]; vDSP_DFT_ExecuteD(inverse, real, imag, hilbertreal, hilbertimag); // Now we have the completed hilbert transform up to a scale factor of n. // We can unscale using vDSP_vsmulD. double scale = 1.0/n; vDSP_vsmulD(hilbert, 1, &scale, hilbert, 1, 2*n); printf("Stage 3:\n"); for (int i=0; i<n; ++i) printf("%f%+fi\n", hilbertreal[i], hilbertimag[i]); vDSP_DFT_DestroySetupD(inverse); vDSP_DFT_DestroySetupD(forward); return 0; }
请注意，如果您已经创build了DFT设置并分配了存储，则计算非常简单; “真正的工作”只是：
vDSP_DFT_ExecuteD(forward, even, odd, real, imag); real[0] *= 0.5; real[n/2] = 0.5*imag[0]; imag[0] = 0.0; vDSP_DFT_ExecuteD(inverse, real, imag, hilbertreal, hilbertimag); double scale = 1.0/n; vDSP_vsmulD(hilbert, 1, &scale, hilbert, 1, 2*n);
示例输出：
Stage 2: 36.000000+0.000000i -8.000000+19.313708i -8.000000+8.000000i -8.000000+3.313708i -4.000000+0.000000i 0.000000+0.000000i 0.000000+0.000000i 0.000000+0.000000i Stage 3: 1.000000+3.828427i 2.000000-1.000000i 3.000000-1.000000i 4.000000-1.828427i 5.000000-1.828427i 6.000000-1.000000i 7.000000-1.000000i 8.000000+3.828427i