阶次分析matlab

%~~~~~~~~~~~~~导入数据~~~~~~~~~~~~~%
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（211）plot（tarray_time_amp）title（‘time dominant振动信号时域图‘）xlabel（‘时间time‘）ylabel（‘幅值amplitude‘）;
subplot（211）plot（tarray_time_amp）title（‘变转速数据‘）xlabel（‘时间time‘）ylabel（‘幅值amplitude‘）;
grid on;
% subplot（212）plot（tpluse）title（‘keyphasor键相脉冲仿真信号时域图‘）xlabel（‘时间time‘）ylabel（‘幅值amplitude‘）;
subplot（212）plot（tpluse）title（‘转速脉冲‘）xlabel（‘时间time‘）ylabel（‘幅值amplitude‘）;
grid on;


% %~~~~~~~~~~~~~~~计算脉冲发生时刻~~~~~~~~~~~~~~~%
% Num_pluse1=1;
Num_pluse1=[];
Threshold=1.5;%设定脉冲阈值
for temp2=1:length（pluse）-1;%设temp2为步长为1的[12.....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_pluse1temp2+1];%Num_pluse1=[121967.....temp（2+1）]
    end
end
 Num_pluse1=Num_pluse1（1:length（Num_pluse1））;
if length（Num_pluse1）<2returnend;%如果脉冲信号就一个数则跳出计算脉冲发生时刻的函数
 t_pluse=（Num_pluse1-1）/Fs;%
% 
% 
% 
% 
%~~~~~~~~~~~~~~~等角度时间计算~~~~~~~~~~~~~~~~%
delt_thet=pi/24;
t_angle=[];
for temp3=3:length（t_pluse）;
b=inv（[1t_pluse（temp3-2）t_pluse（temp3-2）^2;1t_pluse（temp3-1）t_pluse（temp3-1）^2;1t_pluse（temp3）t_pluse（temp3）^2]）*[02*pi4*pi]‘;
    if temp3==3;%（如果数temp3等于3则执行第一个循环，如果temp3不等于3则执行第二个循环，即在脉冲时刻[tt0tt3]执行第一个循环，等角度时刻记到[t1，t72]）
       k=0;
           while k<1.5*2*pi/delt_thet;%0                  tt=（sqrt（4*b（3）*（k*delt_thet-b（1））+b（2）^2）-b（2））/（2*（b（3）+eps））;
                  t_angle=[t_anglett];
                  k=k+1;
            end
    else
        k=pi/delt_thet;%k=24（如果temp3不等于3则执行第二个循环，即在脉冲时刻[tt1tt4].......[tt13tt16]执行第二个循环，等角度时刻记到[t72，t120]，[t120t168]........[t648t696]））
            while k>=pi/delt_thet && k<1.5*2*pi/delt_thet;% 24                tt=（sqrt（4*b（3）*（k*delt_thet-b（1））+b（2）^2）-b（2））/（2*（b（3）+eps））;
                t_angle=[t_anglett];
                k=k+1;
            end
    end
end
ar