用希爾伯特黃變換(HHT)求時(shí)頻譜和邊際譜

1、精品文檔用希爾伯特黃變換(hht求時(shí)頻譜和邊際譜1 .什么是hht？hht1是先將信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解(emd解)，然后將分解后的每個(gè)imf分量進(jìn)行hilbert變換，得到信號的時(shí)頻屬性的一種時(shí)頻分析方法。2 .emd分解的步驟。“防解的流程圖如下：3 . 實(shí)例演示。給定頻率分別為 10hz 和 35hz 的兩個(gè)正弦信號相疊加的復(fù)合信號， 采樣頻率 fs=2048hz 的信號，表達(dá)式如下： y=5sin(2*pi*10t)+5*sin(2*pi*35t)(1) 為了對比，先用 fft 對求上述信號的幅頻和相頻曲線。function fftfenxiclear;clc;n=2048;%fft 默

2、認(rèn)計(jì)算的信號是從0 開始的t=linspace(1,2,n);deta=t(2)-t(1);1/detax=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&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)%

3、將計(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=ta3n=2048;%fft 默認(rèn)計(jì)算的信號是從0 開始的t=lin

4、space(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.a2+xi.a2);%計(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ù)制代碼hilbert 變換是求?。ǔ霈F(xiàn)負(fù)頻，實(shí)際小結(jié)：傅里葉變換不能得到瞬時(shí)頻率，即不能得到某個(gè)時(shí)刻的頻率值。瞬時(shí)頻率的方法，但如

5、果只用 hilbert 變換求出來的瞬時(shí)頻率也不準(zhǔn)確。上負(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); endx

6、g;for i=1:m-1%計(jì)算各imf 的方差貢獻(xiàn)率%定義：方差為平方的均值減去均值的平方%均值的平方%imfp2=mean(c(i,:),2).a2%平方的均值%imf2p=mean(c(i,:).a2,2)%各個(gè)imf 的方差mse(i)=mean(c(i,:).a2,2)-mean(c(i,:),2).a2;end;mmse=sum(mse);for i=1:m-1mse(i)=mean(c(i,:).a2,2)-mean(c(i,:),2).a2;%方差百分比，也就是方差貢獻(xiàn)率mseb(i)=mse(i)/mmse*100;%顯示各個(gè)imf 的方差和貢獻(xiàn)率end;%畫出每個(gè)imf 分

7、量及最后一個(gè)剩余分量residual 的圖形figure(1)for i=1:m-1disp(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-1subplot(m+1,1,i+1);set(gcf,color,w)plot(t,c(i,:),k)set(gca,fontname,times new roman)set(gca,fontsi

8、ze,14.0)ylabel(imf,int2str(i)endsubplot(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,

9、fontsize,14.0)ylabel(initial signal,int2str(m-1),amplitude)for i=1:m-1subplot(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)endsubplot(m+1,1,m+1);set(gcf,color,w)f,z=fftfenxi(t,c(m,:);plot(f,z,k)set(

10、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.a2+xi.a2);%計(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í)頻率 )%十算hht寸頻譜和邊際譜a,fa,tt=hhspectrum(c);e,tt1=toimage(a,f

11、a,tt,length(tt);figure(3)disp_hhs(e,tt1) % 二維圖顯示hht時(shí)頻譜，e是求得的hht譜 pausefigure(4)for i=1:size(c,1)faa=fa(i,:);fa,tt1=meshgrid(faa,tt1);%三維圖顯示 hht寸頻圖surf(fa,tt1,e)title(hht 時(shí)頻譜三維顯示)hold onendhold offe=flipud(e);for k=1:size(e,1)bjp(k)=sum(e(k,:)*1/fs;endf=(1:n-2)/n*(fs/2);figure(5)plot(f,bjp);xlabel( 頻

12、率 / hz);ylabel( 信號幅值 );title( 信號邊際譜)要求邊際譜必須先對信號進(jìn)行em防解function a,f,tt = hhspectrum(x,t,l,aff)error(nargchk(1,4,nargin);if nargin 2t=1:size(x,2);endif nargin 3l=1;endif nargin 4aff = 0;endif min(size(x) = 1if size(x,2) = 1x = x;if nargin = 0error(inf doit etre 0)endm=max(max(im);im = log10(im/m+1e-300

13、);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=2anextpow2(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

14、就是 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)，就表示一定有該頻率的振動波出現(xiàn)， 也就是說， 邊際譜能比較準(zhǔn)確地反映信號的實(shí)際頻率成分。而傅里葉變換只能處理平穩(wěn)信號。(2) hht 與 hilbert