C ++中的Matlab希尔伯特变换答案

【问题标题】：Matlab Hilbert Transform in C++C ++中的Matlab希尔伯特变换
【发布时间】：2012-08-10 04:29:28
【问题描述】：

首先，请原谅我在这个领域的无知，我是一名程序员，但一直陷入超出我专业知识（数学和信号处理）的境地。

我有一个需要移植到 C++ 程序的 Matlab 脚本（无需将 matlab 代码编译成 DLL）。它使用带有一个参数的hilbert() 函数。我正在尝试找到一种在 C++ 中实现相同功能的方法（即，有一个函数也只接受一个参数，并返回相同的值）。

我已经阅读了使用 FFT 和 IFFT 构建它的方法，但似乎无法获得像 Matlab 版本那样简单的任何东西。主要的是我需要它来处理 128*2000 矩阵，而我在搜索中发现的任何内容都没有告诉我如何做到这一点。

我可以返回一个复杂的值，或者只返回一个绝对值。集成到代码中越简单越好。

谢谢。

【问题讨论】：

虽然技术上是一个编程问题，但这个问题更适合DSP 站点。

标签： c++ matlab signal-processing

【解决方案1】：

MatLab 函数 hilbert() 实际上并不直接计算希尔伯特变换，而是计算解析信号，这是大多数情况下需要的。它通过采用 FFT、删除负频率（将阵列的上半部分设置为零）并应用逆 FFT 来实现。如果您有一个不错的 FFT 实现，那么在 C/C++（三行代码）中将是直截了当的。

【讨论】：

是的，我在 Matlab 页面 (mathworks.com/help/toolbox/signal/ref/hilbert.html) 上找到了该算法。我安装了 FFTW 库 (fftw.org)，然后按照算法步骤进行操作。 FFTW 在逆变换上的缩放让我有点困惑，但现在一切正常。
我让它工作得很好，但问题是它很慢。我的意思是，非常慢，比 MATLAB 慢。我正在使用 FFTW 的高级界面以相同的计划执行许多转换，但每次执行 2078 个数据点的 9088 次希尔伯特转换（三个步骤）需要 20 秒。 MATLAB，在同一台计算机上，需要 6 秒。任何想法如何更快地做到这一点？
FFT（快速傅里叶变换）将数据集分解为更小的集合（质数长度）并对每个集合执行 DFT。只有当它可以尽可能多地分解数据时，它才会很快，理想情况下是 2 的幂（最小的素数）。您需要将数据补零到可以轻松分解的长度（例如 4096）。您的数字 2078 会分解为 2 * 1039，这是导致性能不佳的原因。
谢谢，这很有帮助。我将样本数量更改为 2048，并实现了 40% 的时间减少。不幸的是，它仍然比 Matlab 慢得多。还有其他很酷的提示吗？
另外，仅供参考，我发现 2187 以 3 为基数（只有 109 个填充零）比以 2 为基数 4096（2018 填充零）要快得多.

【解决方案2】：

This 看起来还不错，只要你能处理 GPL 许可。 much larger numerical computing resource 的一部分。

【讨论】：

我看过那个代码，不能推荐它。它通过直接卷积计算希尔伯特变换。这种方法有两个问题。首先，与基于 FFT 的快速卷积相比，直接卷积非常慢，特别是因为希尔伯特变换的核衰减非常慢。其次，HT核在零处有一个奇异点，卷积积分的数值计算是病态的。
@AndréBergner 感谢您的意见。但至于速度，执行频域卷积所需的 FFT、IFFT 所需的时间呢？这样总共的时间不是会更大吗？
不，这种技术称为快速卷积，在音频处理中非常流行，例如实时对音乐应用长混响。而且，正如我所写，奇点还存在数字问题。在那个链接的软件中，他们通过稍微移动内核来避免它，这样奇点就在两个样本之间。但这会给输出增加分数延迟，这可能是不希望的，因为它在构建分析信号时会扭曲结果。顺便说一句，还有另一种基于快速递归 (IIR) 滤波器的快速逼近方法——所谓的希尔伯特变换器。

【解决方案3】：

下面的简单代码。（注意：这是一个更大项目的一部分）。 L 的值基于您对订单 N 的确定。其中 N = 2L-1。将 N 舍入为奇数。下面的 xbar 基于您定义为设计系统输入的信号。这是在 MATLAB 中实现的。

L = 40;
n = -L:L; % index n from [-40,-39,....,-1,0,1,...,39,40];
h = (1 - (-1).^n)./(pi*n); %impulse response of Hilbert Transform
h(41) = 0; %Corresponds to the 0/0 term (for 41st term, 0, in n vector above)

xhat = conv(h,xbar); %resultant from Hilbert Transform H(w);

plot(abs(xhat))

【讨论】：

【解决方案4】：

这不是您问题的真正答案，但可能是一种让您睡得更好的方法。我相信在矩阵上基本上是 ffts 的特殊情况下，您将无法比 Matlab 快得多。这就是 Matlab 擅长的地方！

Matlab FFT 是使用 FFTW 计算的，FFTW 是用 C 语言编写的事实上最快的 FFT 算法，似乎也由 Matlab 并行化。最重要的是，引用http://www.mathworks.com/help/matlab/ref/fftw.html：

对于 2 的幂的 FFT 维度，介于 214 和 222 之间，MATLAB 软件在其内部数据库中使用特殊的预加载信息优化 FFT 计算。

所以如果你的代码稍微慢一点也不要难过...

【讨论】：