pylb.m

内燃机转子仿真

%实现频域滤波功能
clear all
clc
close all
load 'bneiranji';
load 'bneiranji1';
load 'bneiranji2';
load 'bneiranji3';
load 'bneiranji4';
step=h;
it=7;
ip=2;
if ip==1
y1=aaii(it,:);
y2=aaii1(it,:);
y3=aaii2(it,:);
y4=aaii3(it,:);
y5=aaii4(it,:);
elseif ip==2
y1=xyz(it,:);
y2=xyz1(it,:);
y3=xyz2(it,:);
y4=xyz3(it,:);
y5=xyz4(it,:);
elseif ip==3;
    y1=vv(it,:);
y2=vv1(it,:);
y3=vv2(it,:);
y4=vv3(it,:);
y5=vv4(it,:);
end


ng=length(y1);
nfft=2^nextpow2(ng);

yy1=fft(y1,nfft);
yy2=fft(y2,nfft);
yy3=fft(y3,nfft);
yy4=fft(y4,nfft);
yy5=fft(y5,nfft);

cf=1/step;
t=(0:step:(ng-1)/cf);

fmin=306;
fmax=600;
n1=round(fmin*nfft/cf+1);
n2=round(fmax*nfft/cf+1);
aa=zeros(5,nfft);
aa(1,n1:n2)=yy1(n1:n2);
aa(1,nfft-n2+1:nfft-n1+1)=yy1(nfft-n2+1:nfft-n1+1);
aa(2,n1:n2)=yy2(n1:n2);
aa(2,nfft-n2+1:nfft-n1+1)=yy2(nfft-n2+1:nfft-n1+1);
aa(3,n1:n2)=yy3(n1:n2);
aa(3,nfft-n2+1:nfft-n1+1)=yy3(nfft-n2+1:nfft-n1+1);
aa(4,n1:n2)=yy4(n1:n2);
aa(4,nfft-n2+1:nfft-n1+1)=yy4(nfft-n2+1:nfft-n1+1);
aa(5,n1:n2)=yy5(n1:n2);
aa(5,nfft-n2+1:nfft-n1+1)=yy5(nfft-n2+1:nfft-n1+1);

iy1=ifft(aa(1,:),nfft);
iy2=ifft(aa(2,:),nfft);
iy3=ifft(aa(3,:),nfft);
iy4=ifft(aa(4,:),nfft);
iy5=ifft(aa(5,:),nfft);

iy1=(real(iy1))';
iy2=(real(iy2))';
iy3=(real(iy3))';
iy4=(real(iy4))';
iy5=(real(iy5))';
figure(1);
subplot(6,1,1);
plot(t,y1(1:ng));
subplot(6,1,2);
plot(t,iy1(1:ng));
subplot(6,1,3);
plot(t,iy2(1:ng));
subplot(6,1,4);
plot(t,iy3(1:ng));
subplot(6,1,5);
plot(t,iy4(1:ng));


subplot(6,1,6);
plot(t,iy5(1:ng));
xlabel('sx');
ylabel('sy');
grid on;

figure(2);
subplot(2,1,1);
plot(t,iy1(1:ng));
subplot(2,1,2);
plot(t,iy2(1:ng));

gly1=fft(iy1,nfft);
pgly1=gly1.*conj(gly1)/nfft;

gly2=fft(iy2,nfft);
pgly2=gly2.*conj(gly2)/nfft;

gly3=fft(iy3,nfft);
pgly3=gly3.*conj(gly3)/nfft;

gly4=fft(iy4,nfft);
pgly4=gly4.*conj(gly4)/nfft;

gly5=fft(iy5,nfft);
pgly5=gly5.*conj(gly5)/nfft;

ff=(0:cf/nfft:cf/2-cf/nfft);
nn=1:length(ff);
nn1=length(ff);
subplot(2,1,1);
plot(ff(nn),pgly1(nn));
xlim([fmin,fmax]);
subplot(2,1,2);
plot(ff(nn),pgly2(nn));
xlim([fmin,fmax]);
xxx=ones(1,nn1);
figure(3);
plot3(xxx,ff(nn),pgly1(nn),2*xxx,ff(nn),pgly2(nn),3*xxx,ff(nn),pgly3(nn),4*xxx,ff(nn),pgly4(nn),...
    5*xxx,ff(nn),pgly5(nn));
xlim([1,5]);
ylim([fmin,fmax]);
title('不同裂纹程度轴段5相对扭角功率谱图');
ylabel('频率(Hz)');
set(gca,'XTickLabel',{'无裂纹';'微裂纹';'轻度裂纹';'中度裂纹';'重度裂纹'});