希尔伯特变换（matlab）

1.什么是Hilbert Transform

百度百科
维基百科 （需要梯子科学上网）

2.Analytic Signal of a Sequence 解析序列

clear; clc; close all;
% 生成一个序列，用希尔伯特分析信号
xr = [1 2 3 4];
x = hilbert(xr)

% 其虚部
x_im = imag(x)
% 其实部
x_re = real(x)

% 快速傅里叶变换
x_dft = fft(x)
xr_dft = fft(xr)

x =

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


x_im =

     1    -1    -1     1


x_re =

     1     2     3     4


x_dft =

  10.0000 + 0.0000i  -4.0000 + 4.0000i  -2.0000 + 0.0000i   0.0000 + 0.0000i


xr_dft =

  10.0000 + 0.0000i  -2.0000 + 2.0000i  -2.0000 + 0.0000i  -2.0000 - 2.0000i

3.Analytic Signal and Hilbert Transform 解析信号和希尔伯特变换

% 希尔伯特函数为有限长数据找到精确的解析信号
% 也可以通过使用FIR希尔伯特变换滤波器来生成解析信号，以计算虚部的近似
% 生成一组频率为203，721，1001Hz的正弦序列
% 采样频率为10kHz，时长为1秒
clear; clc; close all;
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)')

% 计算原始序列和解析信号的Welch估计功率谱密度
% 窗长为256，无重叠部分的汉明窗
% 验证解析信号在负频率处没有功率
pwelch([x;y].',256,0,[],fs,'centered')  % help pwelch
legend('Original','hilbert')

% 设计一个60阶希尔伯特变换FIR滤波器，过渡带宽为400Hz
fo = 60;

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

freqz(d,1024,fs)

% 对正弦序列进行滤波以逼近解析信号的虚部
hb = filter(d,x)

% FIR Filter
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)，并与希尔伯特结果进行比较
pwelch([y;[y2 zeros(1,grd)]].',256,0,[],fs,'centered')
legend('hilbert','FIR Filter')

上图为解析信号x(t)的实部与虚部；
计算原始序列和解析信号的Welch估计功率谱密度；解析信号频率大于0
设计一个60阶希尔伯特变换FIR滤波器，过渡带宽为400Hz；
对比前面的图，对正弦序列进行滤波也可以逼近解析信号的虚部；
估计近似解析信号的功率谱密度(PSD)，并与希尔伯特结果进行比较，与前图效果相差不大。