多种算法实现电机故障诊断

clear all clc fs=200; fm=50; u=5; N=400; n=(0:N-1)-N/2; T=N/fs; x=n/fs; f=(-N/2:N/2-1)*fs/N; yz=15*exp(j*2*pi*fm*x+j*pi*3*(x).^2)+1.1*exp(j*2*pi*45*x+j*pi*3*(x).^2)+1.6*exp(j*2*pi*55*x+j*pi*3*(x).^2); %% yzfft=fft(yz); absyzfft=abs(yzfft); absyzfft=fliplr(absyzfft); figure(1); plot(f,absyzfft,'k'); xlabel('频率f/Hz'); ylabel('极大模值/A'); %==========================================================================原函数傅里叶变化的图像 P=0:0.005:4;% fractional power for rf=1:length(P); gamma=P(rf); fy=frft(yz.',gamma); [a1(rf),count(rf)]=max((abs(fy)).^2);%获取不同P值下的分数阶Fourier变换模极大值及其对应下标 end figure(2) plot(P,a1,'k'); xlabel('P值'); ylabel('极大模值/A'); %================================================================Search，获取最大值及对应的P值 b=a1(1); index=1; for i=2:length(a1) if a1(i)>b index=i; b=a1(i); end end maxa1=a1(index); a1(index)=0; for i=1:length(P); if a1(i)>0.5*maxa1 anotherindex=i; end end yz1frft=abs(frft(yz.',P(index))) ; yz2frft=abs(frft(yz.',P(anotherindex))) ; if yz1frft>yz2frft index=index; else index=anotherindex; end k1=P(index);%分数阶因子 yzfrft=frft(yz.',k1); Fx1=fliplr(abs(yzfrft).'); figure(3); plot(f,Fx1,'k'); xlabel('频率f/Hz'); ylabel('极大模值/A'); legend(['变换值P=',num2str(k1)]); %===============================================================分数阶变换后的图形 M=length(Fx1'); Fx1=Fx1(1:M); d1=zeros(M,20);d2=zeros(20,3); l=1; [y3,exf,exa,maxc,y2]=Smaxp(Fx1'); while exa>5 d2(l,1:3)=[exf exa maxc]; [y3,exf,exa,maxc,y2]=Smaxp(y3); l=l+1; end L=l-1; d3=d2; exfm=d2(1,1);H=0; for i=1:L exf=d2(i,1); if abs(exf-exfm)<8 d2(i,1)=0;d2(i,2)=0;d2(i,3)=0; end if d2(i,1)~=0 H=H+1; end end d4=zeros(H,3);h=1; for i=1:L if d2(i,1)~=0 s=d2(i,:); d4(h,:)=s;h=h+1; end end %y5=abs(y3); y5=zeros(M,1); for i=1:H exf=d4(i,1); exa=d4(i,2); maxc=d4(i,3); yr=Ssinp(exa,exf,maxc,M); %yr=abs(yr); y5=y5+yr; end %figure(4),plot(y5) y5=abs(y5); %y4=10*log(y4/300)/log(10); S=400; Fxx=y5(1:S); figure(4);plot(f,Fxx,'k'); xlabel('频率f/Hz'); ylabel('极大模值/A');