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 運作結果
