希尔伯特变换(分析信号)使用Apple的加速框架?
Hilbert Transform (Analytical Signal) using Apple's Accelerate Framework?
我在使用苹果的Accelerate Framework
在C++
中获得等效于Matlab
的Hilbert变换时遇到问题。我已经能够使vDSP的FFT算法发挥作用,并且在Paul R的帖子的帮助下,成功地获得了与Matlab相同的结果。
我已经阅读了这两个问题:Jordan的这个stackoverflow问题,并阅读了Matlab算法(在"算法"子标题下)。将算法总结为三个阶段:
- 对输入进行FFT
- 零反射频率和直流和奈奎斯特之间的双频
- 对修改后的正向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++示例(使用苹果的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"第三阶段"的结果(或者至少是想象的部分),但我不确定如何进行,我已经尝试了我能想到的一切,但都没有成功。我有一种感觉,这是因为我没有将Apple Accelerate版本中的反射频率归零,因为它们没有提供(由于是根据DC和奈奎斯特之间的频率计算的)-所以我认为FFT只是将倍频的共轭值作为反射值,这是错误的。如果有人能帮我克服这个问题,我将不胜感激
代码
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: %.4fn", hilbertResult.at(i*2), hilbertResult.at(i*2+1));
}
上面的代码是我用来获得"Stage 3"C++(使用苹果的Accelerate Framework)结果的代码。用于前向fft的fft(..)
方法执行vDSP fft例程,然后执行0.5的标度,然后解压缩(根据Paul R的帖子)。当执行逆fft时,数据被打包,按2.0缩放,使用vDSP进行fft运算,最后按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%+fin", 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%+fin", hilbertreal[i], hilbertimag[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);
样本输出:
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
- 如何创建一个CMake变量,除非显式重写,否则使用默认值
- C++:TypeDef使用元组
- 使用std::multimap迭代器创建std::list
- 从不同线程使用int64的不同字节安全吗
- 比较并显示使用最小值(a,b)和最大值(a、b)升序排列的4个数字
- 为什么在全局范围内使用"extern int a"似乎不行?
- 在C#中处理C++指针而不使用unsafe的最佳方法
- 使用C++库在Android项目中修改gradle中的cmake参数,用于插入指令的测试
- std::variant<>::get() 不能使用 Apple LLVM 10.0 编译
- Apple Mach-O linker 错误使用 JUCE 和重型编译器
- 使用 Apple 的 LLVM 编译器编译 -O 时C++代码段错误,但不使用 g++ -7.2.0
- Apple Mach-O Linker (ld) 错误.ld: -r 和 -dead-strip 不能一起使用
- 在Apple llvm 4.1上使用libc++与libstc++时,c++11支持的确切区别是什么
- C++Mac上使用Linux库(x86_64-Apple-Darwin上的elf64-x86-64)
- 希尔伯特变换(分析信号)使用Apple的加速框架?
- 使用Apple Accelerate Framework vForce库来提高性能
- Apple LLVM 4.2 段错误使用基于范围的循环和引用
- Apple Swift中使用的C++加密库
- 简单的 std::regex_search() 代码无法使用 Apple clang++ -std=c++14 编译
- 在Visual Studio 2015中使用C++制作的跨平台应用程序中,是否可以从Apple的Appkit使用NSEvent?