%~~~~~~~~~~~~~導入數據~~~~~~~~~~~~~%
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
a