使用苹果的Accelerate框架的希尔伯特变换（分析信号）？ [英] Hilbert Transform (Analytical Signal) using Apple's Accelerate Framework?
问题描述
我在使用Apple的 Matlab 等效Hilbert变换时遇到问题>加速框架。我已经能够得到vDSP的FFT算法工作，并在 Paul R的帖子的帮助下，设法得到了结果为Matlab。
我读过两个：这个 stackoverflow问题由Jordan ，并已阅读 Matlab算法（在算法标题）。要将算法总计为3个阶段：
- 进行输入FFT。
- 零反射
- Matlab：
m = [1 2 3 4 5 6 7 8] - C ++：
float m [] = {1,2,3,4,5,6,7,8}; - Take forward FFT of input.
- Zero reflection frequencies and double frequencies between DC and Nyquist.
- Take inverse FFT of the modified forward FFT output.
- Matlab:
m = [1 2 3 4 5 6 7 8]; - C++:
float m[] = {1,2,3,4,5,6,7,8};
下面是每个阶段的Matlab和C ++的输出。示例使用以下数组：
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.6569 i
-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框架）
/ em>
实际：36.0000，Imag：0.0000
实际：-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 ：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的Stage 3结果（或至少虚部），但我不确定如何去做，我试过了我可以想到的一切，没有成功。我有一种感觉，因为我不会清零反射频率在苹果加速版本，因为他们没有提供（由于从直流和奈奎斯特之间的频率计算） - 所以我认为FFT只是采取的共轭的双倍频率作为反射值，这是错误的。如果任何人可以帮助我克服这个问题，我会非常感谢它。
代码
void hilbert（std :: vector< float>& data，std :: vector< float>& hilbertResult ）{
//设置变量。
dataSize_2 = data.size（）* 0.5;
//清除和调整向量。
workspace.clear（）;
hilbertResult.clear（）;
workspace.resize（data.size（）/ 2 + 1，0.0f）;
hilbertResult.resize（data.size（），0.0f）;
//将数据复制到hilbertResult向量中。
std :: copy（data.begin（），data.end（），hilbertResult.begin（））;
//执行正向FFT。
fft（hilbertResult，hilbertResult.size（），FFT_FORWARD）;
//用1s和2s填充工作区（DC和Nyquist之间的双频率）。
workspace.at（0）= workspace.at（dataSize_2）= 1.0f;
for（i = 1; i
//通过工作区向量乘法转发FFT输出。
for（i = 0; i
hilbertResult.at（i * 2 + 1）* = workspace.at（i）;
}
//执行逆FFT。
fft（hilbertResult，hilbertResult.size（），FFT_INVERSE）;
for（i = 0; i
hilbertResult.at（i * 2），hilbertResult.at（i * 2 + 1））;
}
上面的代码是我用来获取'Stage 3'C ++与苹果的加速框架）结果。用于forward fft的 fft（..）方法执行vDSP fft例程，然后是0.5的比例，然后解包（根据Paul R的帖子）。当执行逆fft时，数据被打包，使用vDSP缩放2.0，fft'd，最后缩放1 /（2 * N）。
因此，主要问题（据我所知，因为您的代码示例不包含对vDSP的实际调用）是您尝试使用实数到复杂的FFT例程对于步骤三，这基本上是一个复数到复数的反向FFT（由您的Matlab结果具有非零的虚部）。
这是一个简单C实现使用vDSP匹配您期望的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）;
//看起来像错字？实数到复数DFT将其输入分为
//偶数和奇数索引元素。由于实信号是[1，2，3，...]，
//信号[0]是1，信号[2]是3，对偶数索引等等。
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）;
//在这一点上，我们有前向实数到复数DFT，它与
// MATLAB一致，最多为2。因为除了DC和NY
//之外，除了Hilbert变换的一部分外，我们还需要将所有的值都加倍，所以我不打算对
//进行换算，值
//我们真的想要的。因此，我们只需要将NY移动到正确的地方，
//并将DC和NY缩放0.5。反射频率已经是
//归零，因为实数到复数DFT只写入真实和imag的第一个n / 2
//元素。
real [0] * = 0.5; real [n / 2] = 0.5 * imag [0]; imag [0] = 0.0;
printf（Stage 2：\\\
）;
for（int i = 0; i
double hilbert [2 * n];
double * hilbertreal =& hilbert [0];
double * hilbertimag =& hilbert [n];
vDSP_DFT_ExecuteD（inverse，real，imag，hilbertreal，hilbertimag）;
//现在我们已经完成了希尔伯特变换，直到缩放因子为n。
//我们可以使用vDSP_vsmulD解除缩放。
double scale = 1.0 / n; vDSP_vsmulD（hilbert，1，& scale，hilbert，1,2 * n）;
printf（Stage 3：\\\
）;
for（int i = 0; i
vDSP_DFT_DestroySetupD（inverse）;
vDSP_DFT_DestroySetupD（forward）;
return 0;请注意，如果您已经构建了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）;
输出示例：
第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
阶段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.000000 i
8.000000 + 3.828427i
I am having issues with getting a Matlab equivalent Hilbert transform in C++ with using Apple's Accelerate Framework. I have been able to get vDSP's FFT algorithm working and, with the help of Paul R's post, have managed to get the same outcome as Matlab.
I have read both: this stackoverflow question by Jordan and have read the Matlab algorithm (under the 'Algorithms' sub-heading). To sum the algorithm up in 3 stages:
Below are the outputs of both Matlab and C++ for each stage. The examples use the following arrays:
Matlab Example
Stage 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
Stage 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
Stage 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++ Example (with Apple's Accelerate Framework)
Stage 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
Stage 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
Stage 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
It's quite clear that the end results are not the same. I am looking to be able to compute the Matlab 'Stage 3' results (or at least the imaginary parts) but I am unsure how to go about it, I've tried everything I can think of with no success. I have a feeling that it's because I'm not zeroing out the reflection frequencies in the Apple Accelerate version as they are not provided (due to being calculated from the frequencies between DC and Nyquist) - so I think the FFT is just taking the conjugate of the doubled frequencies as the reflection values, which is wrong. If anyone could help me overcome this issue I would greatly appreciate it.
Code
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));
}
The code above is what I used to get the 'Stage 3' C++ (with Apple's Accelerate Framework) result. The fft(..) method for the forward fft performs the vDSP fft routine followed by a scale of 0.5 and then unpacked (as per Paul R's post). When the inverse fft is performed, the data is packed, scaled by 2.0, fft'd using vDSP and finally scaled by 1/(2*N).
So the main problem (as far as I can tell, as your code sample doesn't include the actual calls into vDSP) is that you’re attempting to use the real-to-complex FFT routines for step three, which is fundamentally a complex-to-complex inverse FFT (as evidenced by the fact that your Matlab results have non-zero imaginary parts).
Here’s a simple C implementation using vDSP that matches your expected Matlab output (I used the more modern vDSP_DFT routines, which should generally be preferred to the older fft routines, but otherwise this is very similar to what you’re doing, and illustrates the need for a real-to-complex forward transform, but complex-to-complex inverse transform):
#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;
}
Note that if you already have your DFT setups built and your storage allocated, the computation is extremely straightforward; the "real work" is just:
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);
Sample output:
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
这篇关于使用苹果的Accelerate框架的希尔伯特变换（分析信号）？的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
