資源描述：
《希爾伯特黃變換算例2.doc》由會員上傳分享，免費在線閱讀，更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫。
1、電力工程信號處理應(yīng)用希爾伯特黃變換【目的】1．了解希爾伯特黃變換的理論知識及應(yīng)用領(lǐng)域2．用Matlab軟件仿真，驗證希爾伯特黃變換的優(yōu)點【希爾伯特黃變換】希爾伯特黃變換（Hilbert-Huangtransform,HHT）首先采用EMD方法將信號分解為若干個IMF分量之和，然后對每個IMF分量進(jìn)行Hilbert變換得到的瞬時頻率和瞬時幅值，從而得到信號的Hilbert譜，Hilbert譜表示了信號完整的時間-頻率分布，是具有一定的自適應(yīng)的時頻分析方法。與前面的小波分析方法相比，避免了小波分析基選取的困難。分析非線性、非平穩(wěn)信號采用基于經(jīng)驗?zāi)?/p>
2、態(tài)分解的HHT方法可以較好地分析信號的局域動態(tài)行為和特征。由于HHT方法的種種特點，其在機械振動、生物醫(yī)學(xué)、故障診斷、海洋學(xué)科、地震工程學(xué)以及經(jīng)濟學(xué)各學(xué)科中得到了廣泛應(yīng)用。在電力系統(tǒng)領(lǐng)域中，HHT方法可用于諧波分析、同步電機參數(shù)辨識、低頻震蕩分析、電能質(zhì)量檢測、磁鐵諧振過電壓辨識等方面和超高速方向保護等方面。HHT方法在電力系統(tǒng)中的應(yīng)用還在進(jìn)一步的研究和探索中。9【EMD分解】對于一個時間序列,其經(jīng)驗?zāi)B(tài)分解過程如下：（1）確定原始信號的所有極大值點和極小值點；（2）采用樣條函數(shù)求出的上、下包絡(luò)線，并計算均值；（3）做差；（4）是否滿足終止條
3、件，若不滿足將作為新的輸入信號轉(zhuǎn)至第（1）步，否則轉(zhuǎn)為第（5）步；（5）令，即為一個IMF分量，做差；（6）是否滿足終止條件，若不滿足則將作為新的輸入信號轉(zhuǎn)至第（1）步，若滿足則EMD分解過程結(jié)束，不能提取的為殘余量。具體流程如圖1所示。9圖1EMD分解流程圖對于分解總階數(shù)為的時間序列，最后可以表示成式中，為殘余函數(shù)，它是以單調(diào)函數(shù)?！舅憷?】考察兩個函數(shù)（1）（2）編程：[EMD分解程序]1.emd.mfunctionimf=emd(x)x=transpose(x(:));imf=[];while~ismonotonic(x)x1=x;sd
4、=Inf;while(sd>0.1)
5、
6、~isimf(x1)s1=getspline(x1);s2=-getspline(-x1);x2=x1-(s1+s2)/2;sd=sum((x1-x2).^2)/sum(x1.^2);x1=x2;endimf{end+1}=x1;x=x-x1;endimf{end+1}=x;9functionu=ismonotonic(x)u1=length(findpeaks(x))*length(findpeaks(-x));ifu1>0u=0;elseu=1;endfunctionu=isimf(x)N=leng
7、th(x);u1=sum(x(1:N-1).*x(2:N)<0);u2=length(findpeaks(x))+length(findpeaks(-x));ifabs(u1-u2)>1u=0;elseu=1;endfunctions=getspline(x)N=length(x);p=findpeaks(x);s=spline([0pN+1],[0x(p)0],1:N);1.Findpesks.mfunctionn=findpeaks(x)n=find(diff(diff(x)>0)<0);u=find(x(n+1)>x(n));n(u)=
8、n(u)+1;2.FFTAnalysis.mfunction[Y,f]=FFTAnalysis(y,Ts)Fs=1/Ts;L=length(y);NFFT=2^nextpow2(L);y=y-mean(y);Y=fft(y,NFFT)/L;Y=2*abs(Y(1:NFFT/2+1));f=Fs/2*linspace(0,1,NFFT/2+1);endfunction[yenvelope,yf,yh,yangle]=HilbertAnalysis(y,Ts)yh=hilbert(y);yenvelope=abs(yh);yangle=unwra
9、p(angle(yh));yf=diff(yangle)/2/pi/Ts;9end1.pot_hht.mfunctionplot_hht(x,imf,Ts)%PlottheHHT.%::Syntax%ThearrayxistheinputsignalandTsisthesamplingperiod.%Exampleonuse:[x,Fs]=wavread('Hum.wav');%plot_hht(x(1:6000),1/Fs);%Func:emd%imf=emd(x);fork=1:length(imf)b(k)=sum(imf{k}.*im
10、f{k});th=unwrap(angle(hilbert(imf{k})));%相位d{k}=diff(th)/Ts/(2*pi);%瞬時頻率end[u,v]=s
