%~~~~~~~~~~~~~导入数据~~~~~~~~~~~~~%
clc
clear
Fs=71680;
N=Fs*5;
t=(0:N-1)/Fs;
Adc = 1; %直流分量幅度
% S=sin(2*pi*t.^2- pi/6)+sin(4*pi*t.^2- pi/6)+sin(8*pi*t.^2- pi/6)
% S2=Adc+ sin(2*pi* t.^2- pi/6);%参考轴的转速为n(t)=60t r/min
S=sin(2*pi*t.^2)+sin(4*pi*t.^2)+sin(8*pi*t.^2);
S2=Adc+ sin(2*pi* t.^2);
array_time_amp=S; %导入时域振动信号
pluse=S2; %导入脉冲信号
figure(1);
% subplot(2,1,1),plot(t,array_time_amp),title('time dominant振动信号时域图'),xlabel('时间time'),ylabel('幅值amplitude');
subplot(2,1,1),plot(t,array_time_amp),title('变转速数据'),xlabel('时间time'),ylabel('幅值amplitude');
grid on;
% subplot(2,1,2),plot(t,pluse),title('keyphasor键相脉冲仿真信号时域图'),xlabel('时间time'),ylabel('幅值amplitude');
subplot(2,1,2),plot(t,pluse),title('转速脉冲'),xlabel('时间time'),ylabel('幅值amplitude');
grid on;
% %~~~~~~~~~~~~~~~计算脉冲发生时刻~~~~~~~~~~~~~~~%
% Num_pluse1=1;
Num_pluse1=[];
Threshold=1.5;%设定脉冲阈值
for temp2=1:length(pluse)-1;%设temp2为步长为1的[1,2.....,n-1]n-1维数组;即将所有的采样点编号;length函数功能是返回pluse的数组维数;for循环体循环n-1次,再结束。
% for temp2=1:length(pluse)-1;
if(pluse(temp2)<1.5&&pluse(temp2+1)>=1.5)
Num_pluse1=[ Num_pluse1,temp2+1];%Num_pluse1=[1,219,6,7.....temp(2+1)]
end
end
Num_pluse1=Num_pluse1(1:length(Num_pluse1));
if length(Num_pluse1)<2,return,end;%如果脉冲信号就一个数则跳出计算脉冲发生时刻的函数
t_pluse=(Num_pluse1-1)/Fs;%
%
%
%
%
%~~~~~~~~~~~~~~~等角度时间计算~~~~~~~~~~~~~~~~%
delt_thet=pi/24;
t_angle=[];
for temp3=3:length(t_pluse);
b=inv([1,t_pluse(temp3-2),t_pluse(temp3-2)^2;1,t_pluse(temp3-1),t_pluse(temp3-1)^2;1,t_pluse(temp3),t_pluse(temp3)^2])*[0,2*pi,4*pi]';
if temp3==3;%(如果数temp3等于3则执行第一个循环,如果temp3不等于3则执行第二个循环,即在脉冲时刻[tt0,tt3]执行第一个循环,等角度时刻记到[t1,t72])
k=0;
while k<1.5*2*pi/delt_thet;%0<k<1.5*2*24=72
tt=(sqrt(4*b(3)*(k*delt_thet-b(1))+b(2)^2)-b(2))/(2*(b(3)+eps));
t_angle=[t_angle,tt];
k=k+1;
end
else
k=pi/delt_thet;%k=24,(如果temp3不等于3则执行第二个循环,即在脉冲时刻[tt1,tt4].......[tt13,tt16]执行第二个循环,等角度时刻记到[t72,t120],[t120,t168]........[t648,t696]))
while k>=pi/delt_thet && k<1.5*2*pi/delt_thet;% 24<k<72时,执行循环。否则循环结束。
tt=(sqrt(4*b(3)*(k*delt_thet-b(1))+b(2)^2)-b(2))/(2*(b(3)+eps));
t_angle=[t_angle,tt];
k=k+1;
end
end
end
array_angle=[1:length(t_angle)].*delt_thet;
%
% ~~~~~~~~~~~~~~~~~~~~~~ 数字跟踪滤波器~~~~~~~~~~~~~~~~~~%
i=2;
while i<=length(t_pluse);
ft=2/(t_pluse(i)-t_pluse(i-1));%[t_pluse(i-1),t_pluse(i)]时间段内的采样频率
[b,a]=ellip(6,3,80,ft*25/(Fs/2));%ft*order/(Fs/2),归一化频率,其中1对应Fs/2
array_time_amp(Num_pluse1(i-1):(Num_pluse1(i)-1))=filter(b,a,array_time_amp(Num_pluse1(i-1):( Num_pluse1(i)-1)));%对[Num_pluse1(1),Num_pluse1(2)-1]内振动信号进行数字跟踪低通滤波
i=i+1;
end
%~~~~~~~~~~~~~~~~~~等角度时刻的幅值拟合~~~~~~~~~~~~~~~~%
syms array_angle_amp;
array_angle_amp=[];
array_angle_amp=interp1(t,array_time_amp,t_angle,'spline');%linear%spline%cubic
%~~~~~~~~~~~~评价标准:幅值平方求能量谱~~~~~~~~~~~~~~~%
nlfx=(0:length(t_angle)-1)*48/length(t_angle);
nlpfz=(abs(fft(array_angle_amp))*2*(1/0.5)/length(array_angle_amp)).^2;
pks = findpeaks(nlpfz(1:length(nlpfz)/2),'minpeakheight',0.1,'NPEAK',3)
znl=sum(nlpfz)
figure(5);
plot(nlfx(1:length(nlfx)/2),nlpfz(1:length(nlpfz)/2)),title('time dominant振动信号时域图');
%------------------加窗--------------------------------%
%-------------------------------------------------------%
w=hanning(length(array_angle_amp))';%rectwin 矩形窗 Kaiser 凯瑟窗 Hamming 海明窗 hanning 汉宁窗 bartlett 三角窗 Gausswin 高斯窗 Blackman 布莱克曼窗 Chebwin 切比雪夫窗
array_angle_amp=array_angle_amp.*w;
% ---------------------------------------------------------------%
%-------------- FFT ---------------------------------------%
angle_dom_ffty=abs(fft(array_angle_amp))*2/length(array_angle_amp);
delt_order=2*pi/(length(angle_dom_ffty)*delt_thet);
angle_dom_fx=(0:length(angle_dom_ffty)-1)*delt_order;
FFTy=abs(fft(array_time_amp))*2/length(array_time_amp);
fx=(0:length(array_time_amp)-1)*Fs/length(array_time_amp)/1000;
%---------------2 维图形显示-------------------------------
figure(2);
subplot(2,1,1),plot(fx(1:length(fx)/2),FFTy(1:length(fx)/2)),title('Frequency domain频域');
xlabel('频率f'), ylabel('幅值amplitude');
title('振动仿真信号 幅度-频率曲线图');
subplot(2,1,2),plot(array_angle,array_angle_amp), title('angle dominant振动仿真信号角域图'), xlabel('角度angle /rad'),
ylabel('幅值amplitude');
grid on;
figure(3);
subplot(1,1,1),plot(angle_dom_fx(1:length(angle_dom_fx)/2),angle_dom_ffty(1:length(angle_dom_fx)/2)),
title('order dominant振动仿真信号的阶次域图'), xlabel('阶次order'),ylabel('幅值amplitude');
grid on;
评论6