logo标签列表

在C++中实现Matlab的希尔伯特变换

c++matlabsignal-processing
4

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

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

我已经研究了使用FFT和IFFT构建它的方法，但似乎没有得到与Matlab版本一样简单的东西。主要问题是我需要它能够在一个128*2000的矩阵上运行，而我在搜索中发现的所有内容都没有展示如何做到这一点。

我可以接受复数值或仅绝对值。越简单集成到代码中，就越好。

谢谢。

虽然这是一个技术上的编程问题，但更适合在DSP网站上提问。 - Hristo Iliev
4个回答
9

MatLab函数hilbert()实际上不是直接计算Hilbert变换，而是计算解析信号，在大多数情况下这是需要的。

它通过进行FFT计算，删除负频率（将数组的上半部分设置为零）并应用反FFT来实现。如果您拥有良好的FFT实现，则在C/C++中很容易实现（只需三行代码）。

是的，我在Matlab页面（http://www.mathworks.com/help/toolbox/signal/ref/hilbert.html）上找到了算法。我安装了FFTW库（http://www.fftw.org/），然后按照算法步骤进行操作。FFTW在反向变换中的缩放让我有点困惑，但现在一切都正常了。 - Jordan
我已经把它完美地运行起来了，但问题是它很慢。我的意思是，非常慢，比MATLAB还要慢。我正在使用FFTW的高级接口来执行许多相同计划的变换，但对于每个2078个数据点的9088个希尔伯特变换（三个步骤），需要20秒的时间。在同一台电脑上，MATLAB只需要6秒钟。有什么办法可以更快地完成吗？ - Jordan
1
FFT（快速傅立叶变换）将数据集分解为较小的集合（长度为质数），并对每个集合执行DFT。它只有在尽可能分解数据，理想情况下是2的幂（最小的质数）时才会变得快速。您需要用零填充数据，使其长度易于分解（例如，取4096）。2078可以分解为2 * 1039，这是性能不佳的原因。 - André Bergner
谢谢，这很有帮助。我将样本数量改为2048，时间减少了40%。不幸的是，它仍然比Matlab慢得多。还有其他好的建议吗？ - Jordan
另外，顺便提一下，我发现以3为底数的质数比以2为底数的4096要快得多（只有109个填充零），因为2187。 - Jordan
不错，很有趣。顺便说一下，大多数人认为它只适用于基数2并且必须是这样，但秘密在于小质数。 - André Bergner
2
这看起来很不错，只要你能处理 GPL 许可证。它是一个更大的数值计算资源的一部分。
我看了一下那段代码，不能推荐使用。它通过直接卷积计算希尔伯特变换，这种方法存在两个问题。首先，与基于FFT的快速卷积相比，直接卷积非常慢，特别是因为希尔伯特变换的核衰减非常缓慢。其次，HT核在零点处具有奇异性，卷积积分的数值计算是不适定的。 - André Bergner
@AndréBergner 感谢您的建议。但是就速度而言，进行频域卷积所需的FFT、IFFT所花费的时间如何？这样总体时间不会更长吗？ - Matt Phillips
不，这种技术被称为快速卷积，在音频处理中非常流行，例如实时将长混响应用于音乐。正如我所写的，奇异性也存在数值问题。在那个链接的软件中，他们通过稍微移动核心来避免它，使奇异性位于两个样本之间。但这会向输出添加分数延迟，这可能不是期望的，因为它在构建解析信号时会扭曲结果。顺便说一下，还有另一种基于快速递归（IIR）滤波器的快速近似方法——所谓的希尔伯特变换器。 - André Bergner
0
以下是简单的代码。（注意：这是一个更大项目的一部分）。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))
0

这并不是对你问题的真正回答，但也许可以让你睡得更好。我相信，在基本上是矩阵fft的特定情况下，你不会比Matlab快多少。这就是Matlab的优势所在！

Matlab FFT使用FFTW计算，这是用C编写的事实上最快的FFT算法，似乎也被Matlab并行化了。此外，引用自http://www.mathworks.com/help/matlab/ref/fftw.html：

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

因此，如果你的代码稍微慢一些，也不要感到难过...

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接