【信号处理】Matlab实现希尔伯特-黄变换(图像信号处理)

1 内容介绍

1998年，Norden E. Huang(黄锷:中国台湾海洋学家)等人提出了​​经验模态分解​​方法，并引入了Hilbert谱的概念和Hilbert​​谱分析​​的方法，美国国家航空和宇航局(NASA)将这一方法命名为Hilbert-Huang Transform，简称HHT，即希尔伯特-黄变换。

HHT主要内容包含两部分，第一部分为经验模态分解(Empirical Mode Decomposition，简称EMD)，它是由Huang提出的;第二部分为Hilbert谱分析(Hilbert Spectrum Analysis，简称HSA)。简单说来，HHT处理​​非平稳信号​​的基本过程是:首先利用EMD方法将给定的信号分解为若干固有​​模态​​函数(以Intrinsic Mode Function或IMF表示，也称作本征模态函数)，这些IMF是满足一定条件的分量;然后，对每一个IMF进行Hilbert变换，得到相应的Hilbert谱，即将每个IMF表示在联合的时频域中;最后，汇总所有IMF的Hilbert谱就会得到原始信号的Hilbert谱。

2 部分代码

clear all;

close all

signal=load('Signal.txt');

N=length(signal);

fs=1840; % 采样频率

wp=[50 150]*2/fs; %对原始信号进行滤波处理

ws=[30 200]*2/fs;

rp=1;

rs=15;

[num wn]=cheb2ord(wp,ws,rp,rs);

[b,a]=cheby2(num,rs,wn);

s=filter(b,a,signal);

S=s'; % 将滤波后信号赋值给S

YSS=S; % 滤波信号为始信号

ts=(1/fs); % 采样时间

t=0:ts:ts*(N-1);

%画图－－－滤波后信号

figure(1)

plot(t,s);

set(gcf,'color',[1,1,1]);

grid on;

xlabel('Time(ms)');

ylabel('amplitude(mv)');

title('时域图');

legend('115kHz');

axis([0,0.5,-inf,inf]);

% 设置程序中需要用到的标志位

overf=0; % 是否结束信号分解过程

overimf=0; % IMF的结束标志

%信号处理部分

SD=1; % 设置程序中需要的SD、初值

HH=0.01*ones(1,N); % 设置程序中需要的HH初值

% 进行第一次信号筛选过程

[H,overf,mline]=dob(S,t,ts); % 调用子程序

p=0; % 设置IMF组件的个数初值

% 主程序

% 第一部分：对信号进行IMF处理

condition=1;

while (overf==0) % 如果不满足信号分解结束过程

[sumnum,zeronum]=computernumber(H); % 调用子程序

SD=computerSD(H,HH); % 调用子程序

if ((zeronum==sumnum)|(abs(zeronum-sumnum)==1))&(mline(1:N)==0)

overimf=1;

condition=11;

elseif (0.2<=SD)&(SD<=0.3) % 如果H(n)满足IMF条件，退出筛选过程

overimf=1;

condition=11;

else

condition=condition+1;

end

if condition<=20

overimf=0;

else

overimf=1;

end

% 对IMF进行Hilbert处理

if overimf==1 % 如果分解值满足IMF条件

condition=1; % 重新设置condition的值，以便为下一次筛选做准备

p=p+1; % 累加IMF个数

I(p,:)=H; % 将筛选的IMF赋给I，作为第p个组件，I的每行为一个IMF

SS=YSS; % 将滤波后信号（原始信号）赋给SS

for i=1:p % 用原始信号减去所有筛选出来的IMF组件

II=I(i,:);

SS=SS-II;

end

S=SS; % 原始信号与IMF相减后的余项

% 判断余项是否可以作为最后一项信号分解组件

L=length(S);

m=0; % 以下是获得余项的导数序列

for i=1:L-1

m=m+1;

q(m)=S(i)-S(i+1); % 对余项进行求导

end

% 检查余项是否是单调或常数

if q<0 % 余项是单调递增

overf=1;

elseif q>0 % 余项是单调递减

overf=1;

elseif q==0 % 余项是恒量

overf=1;

else

overf=0;

[H,overf,mline]=dob(S,t,ts); % 不满足分解结束条件，则将余项作为原信号继续进行分解

end

else % 如果分解值不满足IMF条件，则将该分解值作为原始信号重复主程序内的步骤

S=H; % 将本次分解得到的信号作为原始信号

HH=H; % 将本次分解得到的H赋给HH，用来计算SD

[H,overf,mline]=dob(S,t,ts); % 调用子程序

end

end

% 绘制所有IMF组件图

if p==0 % 如果信号不包括任何IMF组件

figure(2)

plot(t,H);

grid on;

xlabel('Time(ms)');

ylabel('Original Signal');

axis([0,0.5,0,inf]);

elseif (p>0) % 如果信号包括 p 个IMF组件

figure(2);

set(gcf,'color',[1 1 1]);

for i=2:p+1

subplot(p+1,1,i-1);

plot(t,I(i-1,:));

grid on;

xlabel('Time(ms)');

ylabel('IMF');

grid on;

end

xlabel('Time');

ylabel('C');

subplot(p+1,1,p);

plot(t,S);

grid on;

xlabel('Time(ms)');

ylabel('Residual');

axis([0,0.5,0,inf]);

end

% 绘制二维/三维时-频图

% 组件为p时

[w,theray,Magy,w_s,energy]=doHilbert(I,ts,t); % 调用子程序 w--原始的频率，w--p经过平滑后的频率

p=length(Magy(:,1)); % Magy的行数

q=length(Magy(1,:)); % Magy的列数

% 以下部分是获取幅值的最大值和最小值

for i=1:p

HMax(i)=Magy(i,1);

HMin(i)=Magy(i,1);

for j=2:q

if Magy(i,j)>HMax(i)

HMax(i)=Magy(i,j);

end

if Magy(i,j)<HMin(i)

HMin(i)=Magy(i,j);

end

end

end

Max=HMax(1);

Min=HMin(1);

for i=2:p

if HMax(i)>=Max

Max=HMax(i);

end

if HMin(i)<=Min

Min=HMin(i);

end

end

% 第i(0<I<=p)个组件的最大能量值和它所对应的时间值

for i=1:p

maxpoint(i)=Magy(i,1);

for j=2:1:q

if (Magy(i,j))>maxpoint(i)

maxpoint(i)=Magy(i,j)^2;

maxpoint_t(i)=t(j);

end

end

end

for i=1:q

Magy_all(i)=sum(Magy(:,i).^2);

end

maxpoint=Magy_all(1);

for i=2:1:q

if (Magy_all(i))>maxpoint

maxpoint=Magy_all(i);

maxpoint_t=t(i);

end

end

maxpoint

maxpoint_t

% 绘制二维hilbert图

% 所有组件的Hilbert图：时间-频率-fuzhi图

figure(3);

set(gcf,'color',[1 1 1]);

e1=subplot(1,1,1); % 设置句柄 e1

for i=p:-1:1

for j=1:q-1

h(i,j)=plot(t(j),w(i,j)/(2*pi)); % 绘制每个点

grid on;

%set(h(i,j),'linestyle','.');

set(h(i,j),'color',(Magy(i,j)-Min)/(Max-Min)*[1 1 1]); % 以黑色为底色，幅值越大越白

set(e1,'color',[0 0 0]); % 设置底色为黑色

hold on

end

end

hold off

xlabel('time(ms)');

ylabel('amplitude');

title('希尔波特谱图');

set(gca, 'color', 'white'); % 去掉坐标轴的背景色

set(gcf, 'color', 'white') ; % 去掉图像的背景色

axis([0 0.5 -0 300]);

% 所有组件的：时间-能量图(瞬时能量谱)

figure(4)

plot(t,Magy_all,'-');

grid on;

set(gcf,'color',[1 1 1]);

xlabel('time(ms)');

ylabel('energy');

title('瞬时能量谱');

axis([0,0.5,0,inf]);

grid on;

% ：频率-能量图（希尔伯特谱）

e_max=energy(1); %求中心频率及其对应的最大能量

for i=2:length(w_s)

if energy(i)>e_max

e_max=energy(i);

f_c=w_s(i)/(2*pi);

end

end

e_max

f_c

figure(5)

plot(w_s/(2*pi),energy);

grid on;

hold on;

set(gcf,'color',[1 1 1]);

xlabel('frequence(kHz)');

ylabel('energy');

title('希尔伯特能量谱');

axis([0,230,0,inf]);

% 绘制三维的hilbert图

figure(6);

set(gcf,'color',[1 1 1]);

for i=1:p

plot3(t(1:N-1),w(i,1:N-1)/(2*pi),Magy(i,1:N-1).^2,'-','color',[1 0 0]*(i/p));

hold on

3 运行结果

4 参考文献

[1]王明阳, 柳征, 周一宇. 基于希尔波特-黄变换的冲击无线电信号检测[J]. 信号处理, 2006, 22(4):4.

博主简介：擅长​​智能优化算法​​、​​神经网络预测​​、​​信号处理​​、​​元胞自动机​​、​​图像处理​​、​​路径规划​​、​​无人机​​、​​雷达通信​​、​​无线传感器​​等多种领域的Matlab仿真，相关matlab代码问题可私信交流。部分理论引用网络文献，若有侵权联系博主删除。