Main Content

hilbert

使用希尔伯特变换的离散时间解析信号

全页折叠

语法

x = hilbert(xr)
x = hilbert(xr,n)

说明

示例

x = hilbert(xr) 从实数数据序列 xr 中返回解析信号 x。如果 xr 是矩阵，则 hilbert 用于求对应于每列的解析信号。

示例

x = hilbert(xr,n) 使用 n 点快速傅里叶变换 (FFT) 计算希尔伯特变换。根据需要，输入数据用零填或截断为长度 n。

示例

全部折叠

打开实时脚本

定义一个序列并使用 hilbert 计算其解析信号。

xr = [1 2 3 4];
x = hilbert(xr)
x = 1×4 complex

   1.0000 + 1.0000i   2.0000 - 1.0000i   3.0000 - 1.0000i   4.0000 + 1.0000i

x 的虚部是 xr 的希尔伯特变换，实部是 xr 本身。

imx = imag(x)
imx = 1×4

     1    -1    -1     1


rex = real(x)
rex = 1×4

     1     2     3     4

x 的离散傅里叶变换 (DFT) 的后半部分为零。（在此示例中，变换的后半部分就是最后一个元素。）fft(x) 的 DC 和奈奎斯特元素是纯实数。

dft = fft(x)
dft = 1×4 complex

  10.0000 + 0.0000i  -4.0000 + 4.0000i  -2.0000 + 0.0000i   0.0000 + 0.0000i
打开实时脚本

hilbert 函数用于求有限数据块的准确解析信号。您也可以使用有限冲激响应 (FIR) 希尔伯特变换滤波器计算虚部的逼近值，从而生成解析信号。

生成一个由频率为 203、721 和 1001 Hz 的三个正弦波组成的序列。该序列的采样频率为 10 kHz，采样时间约为 1 秒。使用 hilbert 函数计算解析信号。对 0.01 秒和 0.03 秒之间的解析信号进行绘图。

fs = 1e4;
t = 0:1/fs:1; 

x = 2.5 + cos(2*pi*203*t) + sin(2*pi*721*t) + cos(2*pi*1001*t);

y = hilbert(x);

plot(t,real(y),t,imag(y))
xlim([0.01 0.03])
legend('real','imaginary')
title('hilbert Function')
xlabel('Time (s)')

计算原始序列和解析信号的功率谱密度的韦尔奇估计值。将序列分成长度为 256 的汉明窗不重叠节。确认在负频率下解析信号没有功率。

pwelch([x;y].',256,0,[],fs,'centered')
legend('Original','hilbert')

使用 designfilt 函数设计一个 60 阶希尔伯特变换器 FIR 滤波器。指定 400 Hz 的过渡带宽度。可视化该滤波器的频率响应。

fo = 60;

d = designfilt('hilbertfir','FilterOrder',fo, ...
       'TransitionWidth',400,'SampleRate',fs); 

freqz(d,1024,fs)

对正弦序列进行滤波以逼近解析信号的虚部。

hb = filter(d,x);

滤波器的群延迟 grd 等于滤波器阶数的一半。对此延迟进行补偿。删除虚部的前 grd 个采样以及实部和时间向量的后 grd 个采样。对 0.01 秒到 0.03 秒之间的结果进行绘图。

grd = fo/2;
   
y2 = x(1:end-grd) + 1j*hb(grd+1:end);
t2 = t(1:end-grd);

plot(t2,real(y2),t2,imag(y2))
xlim([0.01 0.03])
legend('real','imaginary')
title('FIR Filter')
xlabel('Time (s)')

估计逼近解析信号的功率谱密度 (PSD)，并将其与 hilbert 结果进行比较。

pwelch([y;[y2 zeros(1,grd)]].',256,0,[],fs,'centered')
legend('hilbert','FIR Filter')

输入参数

全部折叠

输入信号，指定为实数值向量或矩阵。如果 xr 是复数，则 hilbert 忽略其虚部。

示例: sin(2*pi*(0:15)/16) 指定正弦波的一个周期。

示例: sin(2*pi*(0:15)'./[16 8]) 指定一个双通道正弦信号。

数据类型: single | double

DFT 长度，指定为正整数标量。

数据类型: single | double

输出参量

全部折叠

解析信号，以向量或矩阵形式返回。

详细信息

全部折叠

解析信号

hilbert 从实数据序列中返回一个复螺旋序列，该序列有时称为解析信号。

解析信号 x = xr + jxi 具有实部 xr（原始数据）和虚部 xi（包含希尔伯特变换）。虚部是原始实序列经过 90° 相移后的版本。因此，正弦变换为余弦，反之，余弦变换为正弦。希尔伯特变换后的序列具有与原始序列相同的振幅和频率成分。该变换包括取决于原始序列相位的相位信息。

希尔伯特变换在计算时间序列的瞬时属性（尤其是振幅和频率）时非常有用。瞬时振幅是复希尔伯特变换的振幅；瞬时频率是瞬时相位角的时间变化率。对于纯正弦波，瞬时振幅和频率是恒定的。然而，瞬时相位呈锯齿形状，反映局部相位角在单个周期内是如何线性变化的。对于由多个正弦波混合而成的信号，其属性是在很短的时间范围内（即局部）计算的平均值，该时间范围不超出两到三个数据点。有关示例，请参阅希尔伯特变换与瞬时频率。

参考文献 [1] 描述用于最小相位重构的科尔莫戈洛夫方法，该方法就是对时间序列的频谱密度对数进行希尔伯特变换。工具箱函数 rceps 用于执行这种重构。

算法

序列 xr 的解析信号具有单边傅里叶变换。也就是说，该变换在负频率处消失。为了逼近解析信号，hilbert 计算输入序列的 FFT，将对应于负频率的 FFT 系数替换为零，然后计算结果的逆 FFT。

hilbert 使用一个四步算法：

  1. 计算输入序列的 FFT，将结果存储在向量 x 中。

  2. 创建一个向量 h，其元素 h(i) 具有以下值：

    • 对于 i = 1, (n/2)+1，值为 1

    • 对于 i = 2, 3, … , (n/2)，值为 2

    • 对于 i = (n/2)+2、…、n，值为 0

  3. 计算 x 和 h 的按元素乘积。

  4. 计算在步骤 3 中获得的序列的逆 FFT，并返回结果的前 n 个元素。

此算法最早是在 [2] 中介绍的。该方法假设输入信号 x 是有限数据块。此假设允许函数精确地去除 x 中的频谱冗余。基于 FIR 滤波的方法只能逼近解析信号，但其优点是它们能够连续对数据进行运算。有关使用 FIR 滤波器计算希尔伯特变换的另一个示例，请参阅Single-Sideband Amplitude Modulation。

参考

[1] Claerbout, Jon F. Fundamentals of Geophysical Data Processing with Applications to Petroleum Prospecting. Oxford, UK: Blackwell, 1985.

[2] Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.

[3] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

扩展功能

C/C++ 代码生成
使用 MATLAB® Coder™ 生成 C 代码和 C++ 代码。

GPU 代码生成
使用 GPU Coder™ 为 NVIDIA® GPU 生成 CUDA® 代码。

版本历史记录

在 R2006a 之前推出

全部展开

另请参阅

| |