MATLAB中hilbert函数用法

目录

语法

说明

示例

序列的解析信号

解析信号和希尔伯特变换

        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）分量和奈奎斯特（Nyquist）分量是纯实数。

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赫兹。该序列以10千赫的采样率进行采样，大约为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赫兹的过渡带宽。可视化滤波器的频率响应。

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 — Input signal

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

n — DFT length

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

x — Analytic signal

        解析信号，返回为向量或矩阵。

解析信号

        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时，h(i) = 1
   • 当i = 2, 3, … , (n/2)时，h(i) = 2
   • 当i = (n/2)+2, … , n时，h(i) = 0

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

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

        这个算法最初是在[2]中介绍的。该技术假定输入信号x是一段有限的数据。这个假设允许函数完全去除x中的频谱冗余。基于FIR滤波的方法只能近似解析信号，但它们的优点是可以连续地对数据进行操作。参见单边带幅度调制以获得使用FIR滤波器计算的希尔伯特变换的另一个示例。