【正文】 x1)。xlabel(‘Time/sec’,’fontsize’,8)。ylabel(‘Disp. /um’,’fontsize’,8)。%軸心軌跡figure(2)plot(x1,y1,)。xlabel(’Displacement x/um’,’fontsize’,8)。 ylabel(’ Displacement y/um’,’fontsize’,8)。%對原信號進行重采樣%resampling from 250Hz into 128Hzfs1=1024。dt1=1/fs1。x11=1:N。x2=1:fs/fs1:N。x22=interp1(x11,x1,x2,’spline’)。 Ns=length(x22)。T=dt1*[0:Ns1]。figure(2)plot(T,x22)%小波包分解freq=fs1/(Ns)*(0:Ns/21)。t=wpdec(x22,4,’db44’)。%基于小波包分解的能量比例計算s(1,:)=wprocef(t,[4 0])。p(1,:)=(abs(fft(s(1,:),freq))).^2。s(2,:)=wprocef(t,[4 1])。p(2,:)=(abs(fft(s(2,:),freq))).^2。s(4,:)=wprocef(t,[4 2])。p(4,:)=(abs(fft(s(4,:),freq))).^2。s(3,:)=wprocef(t,[4 3])。p(3,:)=(abs(fft(s(3,:),freq))).^2。s(7,:)=wprocef(t,[4 4])。p(7,:)=(abs(fft(s(7,:),freq))).^2。s(8,:)=wprocef(t,[4 5])。p(8,:)=(abs(fft(s(8,:),freq))).^2。s(6,:)=wprocef(t,[4 6])。p(6,:)=(abs(fft(s(6,:),freq))).^2。s(5,:)=wprocef(t,[4 7])。p(5,:)=(abs(fft(s(5,:),freq))).^2。s(13,:)=wprocef(t,[4 8])。p(13,:)=(abs(fft(s(13,:),freq))).^2。s(14,:)=wprocef(t,[4 9])。p(14,:)=(abs(fft(s(14,:),freq))).^2。s(16,:)=wprocef(t,[4 10])。p(16,:)=(abs(fft(s(16,:),freq))).^2。s(15,:)=wprocef(t,[4 11])。p(15,:)=(abs(fft(s(15,:),freq))).^2。s(11,:)=wprocef(t,[4 12])。p(11,:)=(abs(fft(s(11,:),freq))).^2。s(12,:)=wprocef(t,[4 13])。p(12,:)=(abs(fft(s(12,:),freq))).^2。s(10,:)=wprocef(t,[4 14])。p(10,:)=(abs(fft(s(10,:),freq))).^2。s(9,:)=wprocef(t,[4 15])。p(9,:)=(abs(fft(s(9,:),freq))).^2。ps=sum(p)。pp=p/psttt=[1:16]*32。ttt=[1:16]。bar(ttt,pp,’r’)。ylabel(‘能量比’,’fontsize’, 10, ‘fontweight’, ‘bold’)。xlabel(‘頻帶’, ‘fontsize’, 10, ‘fontweight’, ‘bold’)。 HHT變換用于分析碰摩轉(zhuǎn)子振動信號的程序%%讀數(shù)據(jù)x=load(‘’, ‘r+’)。X=(x/)*1000。x=X(:,2)。y=X(:,3)。dt=。fs=1/dt。N=length(x)。t=dt*[0:N1]。 %時間軸% EMD分解[c,rn,nIMF]=emd(x’)。figure(4)title(‘IMFs and residual’)for i=1。nIMFsubplot(nIMF+1,1,i)plot(t,c(:,i))。axis tight。 ylabel (I,’fontsize’,8)endsubplot(nIMF+1,1,Nimf+1) plot(t,rn)。 axis tight。 ylabel (‘rn’, ‘fontsize’,8)。 xlabel (’Time/sec’,’fontsize’, 8’)。%瞬時頻率譜Fai=[]。w=[]。for i=1:nIMF Hc1=imag(Hilbert(c(:,i)))。 Fai1=atan(Hc1./c(:,i))。 Fai=[Fai Fai1]。 w1=[diff(Fai1)。0]。 w=[w,w1]。endfigure(5)for i=1:nIMF subplot(nIMF,1,i) plot(t,w(:,i)) ylabel (I, ‘fontsize’, 8)。endxlabel (‘Time/sec’, ‘fontsize’, 8)。%Hilbert時頻譜c=c’ 。[A,fa,tt]=hhspectrum(c)。[im,tt]=toimage(A,fa)。disp_hhs(im,[],fs)。%Hilbert邊際譜NN=size(im,1)。 for k=1:NN bjp(k)=sum(im(k,:))*1/fs。 endf=(0:NN1)/NN*(fs/2)。figure(7) plot(f,bjp)。 xlabel(‘Frenquency/Hz’, ‘Fontsize’, 8)。 ylabel(‘Amplitude/um’, ‘Fontsize’, 8)。% %該函數(shù)用來完成EMD分解function[c,rn,Nimf]=emd(y,Ns,fs)s1=y。 % 數(shù)據(jù) Ns=length(y)。 fs=*1000*。s2=s1。for i=1:Nss2(Nsi+1)=s1(i)。ends=[s2 s1 s2]。 %原始信號左右延拓t=0(1/fs):((3*Ns1)/fs)。X0=s。X=X0。tic[Xa Xb]=EvelopeData(X)。tocfigure(3)plot(Xa(:,1),Xa(:,2),’ or’)hold onplot(s)pausepause%EMD分解過程%IMF計算C=[]。h1=X。 IMFflag=1。for i_IMF=1:30 for k=1:420 [Xa Xb]=EvelopData(h1)。 [Ma,Na]=size(Xa)。[Mb,Nb]=size(Xb)。 if min(Ma,Mb)2 Ma Mb IMFflag=0。breakendYa=spline(Xa(:,1)/fs,Xa(:,2),t)。Yb=spline(Xb(:,1)/fs,Xb(:,2),t)。 m1=(Ya+Yb)/2。hik=h1m1。SD=sum(h1h1k).^2/h1.^2)。if SD h1=h1k。 r1=Xh1。 breakendh1=h1k。 end%if IMFflag==0 break。endc1=h1。r1=Xc1。c=[c c1’]。X=r1。h1=X。endc。rn=r1。 %殘余t=0:(1/fs) :((Ns1)/fs)。c=c(Ns+1:2*Ns,:)。rn=rn(Ns+1:2*Ns)。rn=rn’ 。[Ns3,nIMF]=size(c)。nIMF%%該函數(shù)用來計算上下包絡(luò)的最大值和最小值function[Xa,Xb]=ProfileData(X)Ns=length(X)。xmax=0。Xa=[]。xmin=0。Xb=[]。 for i=2:Ns1 if X(i)X(i1) if X(i)X(i+1) Xa=[Xa。[i X(i)]]。 end end endfor i=2:Ns1if X(i)X(i1) if X(i)X(i+1) Xb[Xb。[I X(i)]]。end endend%%該函數(shù)用來計算HilbertHuang spectyumfunction[A,f,tt]=hhspectrum(x)t=1:size(x,2)。l=1。if min(size(x))==1 if size(x,2)==1 x=x’ 。 t=1:size(x,2)。 end Nmodes=1。else Nmodes=size(x,1)。endlt=length(t)。tt=t((l+1):(lt1))。for i=1:Nmodes an(i,:)=Hilbert(x(I,:)’)’ 。 f(i,:)=instfreq(an(i.:)’,tt,1)’。 A=abs(an(:,l+1:end1))。end