资源描述：
《希尔伯特黄变换算例2.doc》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、电力工程信号处理应用希尔伯特黄变换【目的】1．了解希尔伯特黄变换的理论知识及应用领域2．用Matlab软件仿真，验证希尔伯特黄变换的优点【希尔伯特黄变换】希尔伯特黄变换（Hilbert-Huangtransform,HHT）首先采用EMD方法将信号分解为若干个IMF分量之和，然后对每个IMF分量进行Hilbert变换得到的瞬时频率和瞬时幅值，从而得到信号的Hilbert谱，Hilbert谱表示了信号完整的时间-频率分布，是具有一定的自适应的时频分析方法。与前面的小波分析方法相比，避免了小波分析基选取的困难。分析非线性、非平稳信号采用基于经验模
2、态分解的HHT方法可以较好地分析信号的局域动态行为和特征。由于HHT方法的种种特点，其在机械振动、生物医学、故障诊断、海洋学科、地震工程学以及经济学各学科中得到了广泛应用。在电力系统领域中，HHT方法可用于谐波分析、同步电机参数辨识、低频震荡分析、电能质量检测、磁铁谐振过电压辨识等方面和超高速方向保护等方面。HHT方法在电力系统中的应用还在进一步的研究和探索中。9【EMD分解】对于一个时间序列,其经验模态分解过程如下：（1）确定原始信号的所有极大值点和极小值点；（2）采用样条函数求出的上、下包络线，并计算均值；（3）做差；（4）是否满足终止条
3、件，若不满足将作为新的输入信号转至第（1）步，否则转为第（5）步；（5）令，即为一个IMF分量，做差；（6）是否满足终止条件，若不满足则将作为新的输入信号转至第（1）步，若满足则EMD分解过程结束，不能提取的为残余量。具体流程如图1所示。9图1EMD分解流程图对于分解总阶数为的时间序列，最后可以表示成式中，为残余函数，它是以单调函数。【算例2】考察两个函数（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
