1.什么是HHT？
HHT就是先將信號進行經驗模態(tài)分解（EMD分解），然后將分解后的每個IMF分量進行Hilbert變換，得到信號的時頻屬性的一種時頻分析方法。


2.EMD分解的步驟。
?

EMD分解的流程圖如下：
?



3.實例演示。
給定頻率分別為10Hz和35Hz的兩個正弦信號相疊加的復合信號，采樣頻率fs=2048Hz的信號，表達式如下：y=5sin(2*pi*10t)+5*sin(2*pi*35t)

(1)為了對比，先用fft對求上述信號的幅頻和相頻曲線 。

function fftfenxi

clear;clc;

N=2048;

%fft默認計算的信號是從0開始的

t=linspace(1,2,N);deta=t(2)-t(1);1/deta

x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);

% N1=256;N2=512;w1=0.2*2*pi;w2=0.3*2*pi;w3=0.4*2*pi;

% x=(t>=-200&t<=-200+N1*deta).*sin(w1*t)+(t>-200+N1*deta&t<=-200+N2*deta).*sin(w2*t)+(t>-200+N2*deta&t<=200).*sin(w3*t);

y = x;

m=0:N-1;

f=1./(N*deta)*m;%可以查看課本就是這樣定義橫坐標頻率范圍的

%下面計算的Y就是x(t)的傅里葉變換數(shù)值

%Y=exp(i*4*pi*f).*fft(y)%將計算出來的頻譜乘以exp(i*4*pi*f)得到頻移后[-2,2]之間的頻譜值

Y=fft(y);

z=sqrt(Y.*conj(Y));

plot(f(1:100),z(1:100));

title('幅頻曲線')

xiangwei=angle(Y);

figure(2)

plot(f,xiangwei)

title('相頻曲線')

figure(3)

plot(t,y,'r')

%axis([-2,2,0,1.2])

title('原始信號')

復制代碼



（2）用Hilbert變換直接求該信號的瞬時頻率

clear;clc;clf;

%假設待分析的函數(shù)是z=t^3

N=2048;

%fft默認計算的信號是從0開始的

t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;

x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);

z=x;

hx=hilbert(z);

xr=real(hx);xi=imag(hx);

%計算瞬時振幅

sz=sqrt(xr.^2+xi.^2);

%計算瞬時相位

sx=angle(hx);

%計算瞬時頻率

dt=diff(t);

dx=diff(sx);

sp=dx./dt;

plot(t(1:N-1),sp)

title('瞬時頻率')

復制代碼
小結：傅里葉變換不能得到瞬時頻率，即不能得到某個時刻的頻率值。Hilbert變換是求取瞬時頻率的方法，但如果只用Hilbert變換求出來的瞬時頻率也不準確。（出現(xiàn)負頻，實際上負頻沒有意義！）
（3）用HHT求取信號的時頻譜與邊際譜

function HHT

clear;clc;clf;

N=2048;

%fft默認計算的信號是從0開始的

t=linspace(1,2,N);deta=t(2)-t(1);fs=1/deta;

x=5*sin(2*pi*10*t)+5*sin(2*pi*35*t);

z=x;

c=emd(z);

%計算每個IMF分量及最后一個剩余分量residual與原始信號的相關性

[m,n]=size(c);

for i=1:m;

a=corrcoef(c(i,:),z);

xg(i)=a(1,2);

end

xg;

for i=1:m-1

%--------------------------------------------------------------------

%計算各IMF的方差貢獻率

%定義：方差為平方的均值減去均值的平方

%均值的平方

%imfp2=mean(c(i,:),2).^2

%平方的均值

%imf2p=mean(c(i,:).^2,2)

%各個IMF的方差

mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;

end;

mmse=sum(mse);

for i=1:m-1

mse(i)=mean(c(i,:).^2,2)-mean(c(i,:),2).^2;?

%方差百分比，也就是方差貢獻率

mseb(i)=mse(i)/mmse*100;

%顯示各個IMF的方差和貢獻率

end;

%畫出每個IMF分量及最后一個剩余分量residual的圖形

figure(1)

for i=1:m-1

disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);

end;

subplot(m+1,1,1)

plot(t,z)

set(gca,'fontname','times New Roman')

set(gca,'fontsize',14.0)

ylabel(['signal','Amplitude'])

for i=1:m-1

subplot(m+1,1,i+1);

set(gcf,'color','w')

plot(t,c(i,:),'k')

set(gca,'fontname','times New Roman')

set(gca,'fontsize',14.0)

ylabel(['imf',int2str(i)])

end

subplot(m+1,1,m+1);

set(gcf,'color','w')

plot(t,c(m,:),'k')

set(gca,'fontname','times New Roman')

set(gca,'fontsize',14.0)

ylabel(['r',int2str(m-1)])

%畫出每個IMF分量及剩余分量residual的幅頻曲線

figure(2)

subplot(m+1,1,1)

set(gcf,'color','w')

[f,z]=fftfenxi(t,z);

plot(f,z,'k')

set(gca,'fontname','times New Roman')

set(gca,'fontsize',14.0)

ylabel(['initial signal',int2str(m-1),'Amplitude'])

for i=1:m-1

subplot(m+1,1,i+1);

set(gcf,'color','w')

[f,z]=fftfenxi(t,c(i,:));

plot(f,z,'k')

set(gca,'fontname','times New Roman')

set(gca,'fontsize',14.0)

ylabel(['imf',int2str(i),'Amplitude'])

end

subplot(m+1,1,m+1);

set(gcf,'color','w')

[f,z]=fftfenxi(t,c(m,:));

plot(f,z,'k')

set(gca,'fontname','times New Roman')

set(gca,'fontsize',14.0)

ylabel(['r',int2str(m-1),'Amplitude'])

hx=hilbert(z);

xr=real(hx);xi=imag(hx);

%計算瞬時振幅

sz=sqrt(xr.^2+xi.^2);

%計算瞬時相位

sx=angle(hx);

%計算瞬時頻率

dt=diff(t);

dx=diff(sx);

sp=dx./dt;

figure(6)

plot(t(1:N-1),sp)

title('瞬時頻率')

%計算HHT時頻譜和邊際譜

[A,fa,tt]=hhspectrum(c);

[E,tt1]=toimage(A,fa,tt,length(tt));

figure(3)

disp_hhs(E,tt1) %二維圖顯示HHT時頻譜，E是求得的HHT譜

pause

figure(4)

for i=1:size(c,1)

faa=fa(i,:);

[FA,TT1]=meshgrid(faa,tt1);%三維圖顯示HHT時頻圖

surf(FA,TT1,E)

title('HHT時頻譜三維顯示')

hold on

end

hold off

E=flipud(E);

for k=1:size(E,1)

bjp(k)=sum(E(k,:))*1/fs;?

end

f=(1:N-2)/N*(fs/2);

figure(5)

plot(f,bjp);

xlabel('頻率 / Hz');

ylabel('信號幅值');

title('信號邊際譜')%要求邊際譜必須先對信號進行EMD分解

function [A,f,tt] = hhspectrum(x,t,l,aff)

error(nargchk(1,4,nargin));

if nargin < 2

t=1:size(x,2);

end

if nargin < 3

l=1;

end

if nargin < 4

aff = 0;

end

if min(size(x)) == 1

if size(x,2) == 1

x = x';

if nargin < 2

t = 1:size(x,2);

end

end

Nmodes = 1;

else

Nmodes = size(x,1);

end

lt=length(t);

tt=t((l+1):(lt-l));

for i=1:Nmodes

an(i,:)=hilbert(x(i,:)')';

f(i,:)=instfreq(an(i,:)',tt,l)';

A=abs(an(:,l+1:end-l));

if aff

disprog(i,Nmodes,max(Nmodes,100))

end

end

function disp_hhs(im,t,inf)

% DISP_HHS(im,t,inf)

% displays in a new figure the spectrum contained in matrix "im"

% (amplitudes in log).

%

% inputs : - im : image matrix (e.g., output of "toimage")

% - t (optional) : time instants (e.g., output of "toimage")?

% - inf (optional) : -dynamic range in dB (wrt max)

% default : inf = -20

%

% utilisation : disp_hhs(im) ; disp_hhs(im,t) ; disp_hhs(im,inf)?

% disp_hhs(im,t,inf)

figure

colormap(bone)

colormap(1-colormap);

if nargin==1

inf=-20;

t = 1:size(im,2);

end

if nargin == 2

if length(t) == 1

inf = t;

t = 1:size(im,2);

else

inf = -20;

end

end

if inf >= 0

error('inf doit etre < 0')

end

M=max(max(im));

im = log10(im/M+1e-300);

inf=inf/10;

imagesc(t,fliplr((1:size(im,1))/(2*size(im,1))),im,[inf,0]);

set(gca,'YDir','normal')

xlabel(['time'])

ylabel(['normalized frequency'])

title('Hilbert-Huang spectrum')

function [f,z]=fftfenxi(t,y)

L=length(t);N=2^nextpow2(L);

%fft默認計算的信號是從0開始的

t=linspace(t(1),t(L),N);deta=t(2)-t(1);

m=0:N-1;

f=1./(N*deta)*m;

%下面計算的Y就是x(t)的傅里葉變換數(shù)值

%Y=exp(i*4*pi*f).*fft(y)%將計算出來的頻譜乘以exp(i*4*pi*f)得到頻移后[-2,2]之間的頻譜值

Y=fft(y);

z=sqrt(Y.*conj(Y));

復制代碼



4.總結。

（1）邊際譜與傅里葉譜的比較：
? ?? ?? ? 意義不同：邊際譜從統(tǒng)計意義上表征了整組數(shù)據(jù)每個頻率點的累積幅值分布，而傅里葉頻譜的某一點頻率上的幅值表示在整個信號里有一個含有此頻率的三角函數(shù)組分。
? ?? ?? ?作用不同：邊際譜可以處理非平穩(wěn)信號，如果信號中存在某一頻率的能量出現(xiàn)，就表示一定有該頻率的振動波出現(xiàn)，也就是說，邊際譜能比較準確地反映信號的實際頻率成分。而傅里葉變換只能處理平穩(wěn)信號。

（2） HHT與Hilbert變換的比較： ?
? ?? ?? ? Hilbert變換只是單純地求信號的瞬時振幅，頻率和相位，有可能出現(xiàn)沒有意義的負頻率；HHT變換先將信號進行EMD分解，得到的是各個不同尺度的分量，對每一個分量進行Hilbert變換后得到的是有實際意義的瞬時頻率。

總結

以上是生活随笔為你收集整理的hht变换的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。