syjf01.m

内燃机转子仿真

clear all
clc
fin=input('the name of the data:','s');
fid=fopen(fin,'r');
T=fscanf(fid,'%f',1);%计算时间
t1=fscanf(fid,'%f',1);%脉冲力作用时间
t2=fscanf(fid,'%f',1);%驱动力作用时间
n=fscanf(fid,'%d',1);
m=fscanf(fid,'%f',[n,n]);
k=fscanf(fid,'%f',[n,n]);
c=fscanf(fid,'%f',[n,n]);
w=fscanf(fid,'%d',1);
b=fscanf(fid,'%f',1);
F=fscanf(fid,'%d',[n,1]);
xx0=fscanf(fid,'%f',[n,1]);
xx1=fscanf(fid,'%f',[n,1]);
fclose(fid);

stp=60.0/zhuansu/360*2.0;%转动时间步长（2度/秒）
sTd=rem(T,stp);
TT=T-stpd;
nT=TT/stp;%计算步数
st1d=rem(t1,stp)
tt1=T-st1d;
nt1=tt1/stp;%脉冲力作用步数
st2d=rem(t2,stp)
tt2=T-st2d;
nt1=tt2/stp;%驱动力作用步数
%**********************激励计算




N=length(xx(1,:));
for i=1:3
xg(i,2*N-1)=xcorr(xx(i,:));
end
for i=1:3
jz(i,1)=sum(x1(:,i))/N;
jfz(i,1)=sum(x1(:,i).*x1(:,i))/N;
fc(i,1)=sum((x1(:,i)-jz(i,1)).^2)/N;
end
figure(1);
plot(dt*i,xg,'g');title('相关');
figure(2);
plot(dt*i,jz,'g');title('均值');
figure(3);
plot(dt*i,jfz,'g');title('均方值');
figure(4);
plot(dt*i,fc,'g');title('方差');
%频域分析
n=length(xx(:,1));
nfft=2^nextpow2(n);
ft=0:sf/nfft:sf/2-sf/nfft;
for i=1:3
x=(f(:,i))';
y=(y(:,i))'-x;
switch wind
    case 1
    w=boxcar(nfft);
    case 2
     w=hanning(nfft);
    case 3
     w=hamming(nfft);
    case 4
     w=blackman(nfft);
    case 5
     w=triang(nfft);
    otherwise
        w=boxcar(nfft);
end
z1=psd(y,nfft,sf,w,nfft/2);
figure(5);
plot(ft,abs(z1));
xlabel('频率（Hz）');
ylabel('幅值');
z2=csd(x,y,nfft,sf,w,nfft/2);
figure(6);
plot(ft,angle(z2));
xlabel('频率（Hz）');
ylabel('相位');
z3=tfe(x,y,nfft,sf,w,nfft/2);
figure(7);
plot(ft,real(z3));
xlabel('频率（Hz）');
ylabel('实部');
figure(8);
plot(ft,imag(z3));
xlabel('频率（Hz）');
ylabel('虚部');
z4=cohere(x,y,nfft,sf,w,nfft/2);
figure(9);
plot(ft,abs(z4));
xlabel('频率（Hz）');
ylabel('幅值');
end
%时频分析
%倒频谱分析
%高阶谱分析