1.什么是HHT？
HHT就是先將信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解（EMD分解），然后將分解后的每個(gè)IMF分量進(jìn)行Hilbert變換，得到信號的時(shí)頻屬性的一種時(shí)頻分析方法。


2.EMD分解的步驟。
?

EMD分解的流程圖如下：
?



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

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

function fftfenxi

clear;clc;

N=2048;

%fft默認(rèn)計(jì)算的信號是從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;%可以查看課本就是這樣定義橫坐標(biāo)頻率范圍的

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

%Y=exp(i*4*pi*f).*fft(y)%將計(jì)算出來的頻譜乘以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('原始信號')

復(fù)制代碼



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

clear;clc;clf;

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

N=2048;

%fft默認(rèn)計(jì)算的信號是從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);

%計(jì)算瞬時(shí)振幅

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

%計(jì)算瞬時(shí)相位

sx=angle(hx);

%計(jì)算瞬時(shí)頻率

dt=diff(t);

dx=diff(sx);

sp=dx./dt;

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

title('瞬時(shí)頻率')

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

function HHT

clear;clc;clf;

N=2048;

%fft默認(rèn)計(jì)算的信號是從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);

%計(jì)算每個(gè)IMF分量及最后一個(gè)剩余分量residual與原始信號的相關(guān)性

[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

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

%計(jì)算各IMF的方差貢獻(xiàn)率

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

%均值的平方

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

%平方的均值

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

%各個(gè)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;?

%方差百分比，也就是方差貢獻(xiàn)率

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

%顯示各個(gè)IMF的方差和貢獻(xiàn)率

end;

%畫出每個(gè)IMF分量及最后一個(gè)剩余分量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)])

%畫出每個(gè)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);

%計(jì)算瞬時(shí)振幅

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

%計(jì)算瞬時(shí)相位

sx=angle(hx);

%計(jì)算瞬時(shí)頻率

dt=diff(t);

dx=diff(sx);

sp=dx./dt;

figure(6)

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

title('瞬時(shí)頻率')

%計(jì)算HHT時(shí)頻譜和邊際譜

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

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

figure(3)

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

pause

figure(4)

for i=1:size(c,1)

faa=fa(i,:);

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

surf(FA,TT1,E)

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

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('信號邊際譜')%要求邊際譜必須先對信號進(jìn)行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默認(rèn)計(jì)算的信號是從0開始的

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

m=0:N-1;

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

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

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

Y=fft(y);

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

復(fù)制代碼



4.總結(jié)。

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

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。