369 发表于 2014-12-4 10:58

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

本帖最后由 wdhd 于 2016-4-15 15:39 编辑

  %-----数据导入(手动修改)--------------------------------------------
  %--------------------------------------------------------------------------
  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=;
  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));
  =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])*';
  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=;
  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=;
  k=k+1;
  end
  end
  end
  array_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=;
  % x=;
  % angle_amp=interp1(x,y,t_angle(temp5),'linear');
  % array_angle_amp=;
  %
  % % -------- 三次多项式插值拟合------------------
  % % if temp4==1;
  % % y=;
  % % x=;
  % % angle_amp=interp1(x,y,t_angle(temp5),'linear');
  % % array_angle_amp=;
  % % else
  % %
  y=;
  % %由前两点与后两点得用三次多项式插值求得此角度对应的幅值
  % % x=;
  % % angle_amp=interp1(x,y,t_angle(temp5),'cubic');
  % % array_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是什么?
  请前辈们 指导一下 谢谢~~~
页: [1]
查看完整版本: 请教下论坛上阶次跟踪程序的相关问题 谢谢