fs=1024;
N=2048;
t=(0:N-1)./fs;
x1=1.5*(cos(4*pi*(t.^6+(2*t.^3)/3+3*t.^2+30*t))+1);
x2=1.5*cos(10*pi*(t.^6+(2*t.^3)/3+3*t.^2+30*t));
x=x1.*x2;%仿真信号
pluse=t.^5+2*t.^2+6*t+30;%转速信号
Ts = 1/fs;
imf=emd(x);
m=15;
for i=1
A=imf(i,:);
B=abs(A);
PKS=findpeaks(B);
c1=1:1:length(PKS);c2=1:1:length(t);
e=spline(c1,PKS,c2);
f=B./e;
if (all(f(1,:)<1))
break
else
F=f;
end
end
[A,fa,tt]=hhspectrum(F);%求取瞬时频率与振幅:A:是振幅向量,fa:瞬时频率,tt:时间序列号
pluse=fa;%瞬时频率
%%
Num_pluse1=[];Threshold=0.3;
for temp2=1:length(pluse)-1
if(pluse(temp2)<Threshold&&pluse(temp2+1)>=Threshold)
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;
% ~~~~~~~~~~~~~~~~~~~~~~ 数字跟踪滤波器~~~~~~~~~~~~~~~~~~%
i=2;
while i<=length(t_pluse)
f0=2/(t_pluse(i)-t_pluse(i-1));%[t_pluse(i-1),t_pluse(i)]时间段内的采样频率
[b,a]=ellip(60,2,80,f0/(2*fs));%ft*order/(Fs/2),归一化频率,其中1对应Fs/2
c=Num_pluse1(i-1):(Num_pluse1(i)-1);
array_time_amp=filter(b,a,c);%对[Num_pluse1(1),Num_pluse1(2)-1]内振动信号进行数字跟踪低通滤波
i=i+1;
end
%%
%~~~~~~~~~~~~~~~等角度时间计算~~~~~~~~~~~~~~~~%
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<Threshold*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<Threshold*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;
%~~~~~~~~~~~~~~~~~~等角度时刻的幅值拟合~~~~~~~~~~~~~~~~%
% syms array_angle_amp;
% array_angle_amp=[];
array_angle_amp=interp1(t,array_time_amp,t_angle,'cubic');%linear%spline%cubic
x=1:1:10
v=[a1 a2 a3 a4 a5 a6 a7 a8 a8 a10]注释 因为x有十个数,所以v必须有十个数
x1=7.8
v1=interp1(x,v,'spline')