请教下论坛上阶次跟踪程序的相关问题 谢谢

369

　　%-----数据导入(手动修改)--------------------------------------------
　　%--------------------------------------------------------------------------
　　N=2200000;
　　Fs=40000;
　　order_max=48;
　　t=(0:N-1)*1/Fs;
　　array_time_amp=Data1; %时域振动信号导入
　　pluse=Data2; %脉冲信号导入
　　%--------------------------------------------------------------------------
　　%-----计算脉冲发生时刻---------------------------------------------
　　%--------------------------------------------------------------------------
　　Num_pluse1=1;
　　threshold=200; %设定脉冲阈值
　　for temp2=1:length(pluse)-1;
　　if (pluse(temp2)<200&&pluse(temp2+1)>=200)
　　Num_pluse1=[Num_pluse1,temp2+1];
　　end
　　end
　　if length(Num_pluse1)<2, return,end;
　　t_pluse=(Num_pluse1-1)/Fs;
　　%----------------------------------------------------------
　　%-------数字跟踪滤波-----------------------------------
　　%-----------------------------------------------------------
　　i=2;
　　while i<=length(t_pluse);
　　ft=2/(t_pluse(i)-t_pluse(i-1));
　　[b,a]=ellip(6,3,80,ft*order_max/(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)));
　　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_pl
　　use(temp3),t_pluse(temp3)^2])*[0,2*pi,4*pi]';
　　if temp3==3;
　　k=0;
　　while k<1.5*2*pi/delt_thet;
　　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;
　　while k>=pi/delt_thet && k<1.5*2*pi/delt_thet;
　　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;
　　%------等角度时刻的幅值拟合(三种拟合方式)-----
　　%------线性插值拟合-----------------------------
　　% array_angle_amp=[];
　　% for temp5=1:length(t_angle);
　　% for temp4=1:length(t);
　　% if t_angle(temp5)>=t(temp4)&&t_angle(temp5)
　　% %linear interpolation
　　% y=[array_time_amp(temp4),array_time_amp(temp4+1)];
　　% x=[t(temp4),t(temp4+1)];
　　% angle_amp=interp1(x,y,t_angle(temp5),'linear');
　　% array_angle_amp=[array_angle_amp,angle_amp];
　　%
　　% % -------- 三次多项式插值拟合------------------
　　% % if temp4==1;
　　% % y=[array_time_amp(temp4),array_time_amp(temp4+1)];
　　% % x=[t(temp4),t(temp4+1)];
　　% % angle_amp=interp1(x,y,t_angle(temp5),'linear');
　　% % array_angle_amp=[array_angle_amp,angle_amp];
　　% % else
　　% %
　　y=[array_time_amp(temp4-1),array_time_amp(temp4),array_time_amp(temp4+1),array_time_amp(temp4+2)];
　　% %由前两点与后两点得用三次多项式插值求得此角度对应的幅值
　　% % x=[t(temp4-1),t(temp4),t(temp4+1),t(temp4+2)];
　　% % angle_amp=interp1(x,y,t_angle(temp5),'cubic');
　　% % array_angle_amp=[array_angle_amp,angle_amp];
　　% % end
　　% break;
　　% end
　　% end
　　% end
　　%----------三次样条插值拟合-------------------------------
　　array_angle_amp=interp1(t,array_time_amp,t_angle,'cubic');
　　%--------------------------------------------------------
　　%------------------加窗--------------------------------
　　%-------------------------------------------------------
　　w=hanning(length(array_angle_amp))';
　　array_angle_amp=array_angle_amp.*w;
　　s%---------------------------------------------------------------
　　%-------------- 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)/t_sample;
　　%--------------------------------------------------------------
　　%---------------2 维图形显示----------------------------------
　　%----------------------------------------------------------------
　　figure(1);
　　subplot(2,1,1),plot(t,array_time_amp),title('timedomianat'),xlabel('time'),ylabel('amplitude');
　　grid on;
　　subplot(2,1,2),plot(t,pluse),title('keyphasor'),xlabel('time'),ylabel('amplitude');
　　grid on;
　　figure(2);
　　subplot(2,1,1),plot(fx(1:length(fx)/2),FFTy(1:length(fx)/2)),title('Frequencydomain');
　　xlabel('f'), ylabel('amplitude');
　　subplot(2,1,2),plot(array_angle,array_angle_amp), title('angle dominant'),xlabel('angle /rad'),
　　ylabel('amplitude');
　　grid on;
　　figure(3);
　　subplot(2,1,1),plot(angle_dom_fx(1:length(angle_dom_fx)/2),angle_dom_ffty(1:length(angle_dom_fx)/2));
　　1.什么是脉冲阈值?其作用? 运行此程序t_pluse总为0 如何修改?
　　2.重采样过程 维数问题怎么修改?
　　3.fx=(0:length(array_time_amp)-1)/t_sample; t_sample是什么?
　　请前辈们 指导一下 谢谢~~~