Main Content

解析信号和希尔伯特变换

打开实时脚本

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')

另请参阅

|