function sw = checksw
% CHECKSW Check if running MATLAB or GNU/Octave
% sw = checksw
%
% sw 'MATLAB' or 'OCTAVE' depending on what software checksw is
% called from. 0 is returned if neither MATLAB or Octave was
% found (unexpected).
%
% Copyright (c) 2009-2011 by Anders Brandt
% Email: abra@iti.sdu.dk
% Version: 1.0 2011-06-23
% This file is part of ABRAVIBE Toolbox for NVA
v=ver;
sw=0;
for n = 1:length(v)
Name=v(n).Name;
if strcmp(upper(Name),'MATLAB')
sw='MATLAB';
elseif strcmp(upper(Name),'OCTAVE')
sw='OCTAVE';
end
end
function[omega,accfdat,proxfdat,static]=ordtrack(fsamp,keyfsamp,keytdat,acctdat,proxtdat,no_cycles,limits)
% ORDTRACK Order tracking
% [OMEGA,ACCFDAT,PROXFDAT,STATIC]=ORDTRACK(FSAMP,KEYFSAMP,KEYTDAT,...
% ACCTDAT,PROXTDAT,NO_CYCLES,LIMITS);
% returns OMEGA, the running speed in rad/s, ACCFDAT and PROXFDAT,
% the tracked accelerometer and proximitor data and STATIC, the
% static proximitor data given:
% FSAMP, KEYFSAMP - data and keyphasor sampling rates (Hz)
% KEYTDAT - keyphasor time data
% ACCTDAT,PROXTDAT - accelerometer and proximitor time data
% NO_CYCLES - no of cycles over which to average data (default 8). The higher
% the number, the smoother the data (but also the more information
% is lost.)
% LIMITS - [f1,f2], the speeds (in Hz) over which the tracking is
% to be done (default all data);
[pulses,spd,time]=shaftspd(keytdat,keyfsamp);
if ~exist('no_cycles')
no_cycles=[];
end
if isempty(no_cycles)
no_cycles=8;
end
p=fix(pulses/(keyfsamp/fsamp)); % index of pulses in data channels
l=max([size(acctdat,1) size(proxtdat,1)]);
time=(0:l-1)'/fsamp;
if ~exist('limits')
limits=[];
end
if isempty(limits)
x(1)=time(1);
x(2)=time(end);
else
x(1)=limits(1);
x(2)=limits(2);
end
x=sort(x);
time_at_pulses=pulses/fsamp;
[val,x(1)]=min(abs(time_at_pulses-x(1)));
[val,x(2)]=min(abs(time_at_pulses-x(2)));
l1=pulses(x(1));
l2=pulses(x(2));
ll1=fix(l1/(keyfsamp/fsamp));
ll2=fix(l2/(keyfsamp/fsamp));
pp=polyfit(time_at_pulses(x(1):x(2)),spd(x(1):x(2)),2);
fitspd=polyval(pp,time(ll1:ll2));
fitspd_pulses=p(x(1):x(2))-p(x(1))+1;
if ~exist('no_cycles')
disp('Finding pulses for frequency spacing');
[si,ei]=freqspac(fitspd_pulses,fitspd,64,0,res);
else
si=fitspd_pulses(1:no_cycles:length(fitspd_pulses)-no_cycles);
ei=fitspd_pulses(no_cycles+1:no_cycles:length(fitspd_pulses));
end
accfdat=[];
proxfdat=[];
static=[];
if ~isempty(proxtdat)
proxfdat=zeros(length(si),size(proxtdat,2));
for ii=1:size(proxtdat,2)
[freqdat,omega,st]=sinefit(proxtdat(ll1:ll2),time(ll1:ll2),fitspd,si,ei);
proxfdat(:,ii)=freqdat;
static(:,ii)=st;
end
end
if ~isempty(acctdat);
accfdat=zeros(length(si),size(acctdat,2));
for ii=1:size(acctdat,2)
% chno=acc_channels(ii);
[freqdat,omega]=sinefit(acctdat(ll1:ll2),time(ll1:ll2),fitspd,si,ei);
accfdat(:,ii)=freqdat;
end
end
function [rpm,t] = tacho2rpm(x,fs,TLevel,Slope,PPR,NewFs,FilterL)
% TACHO2RPM Extract rpm/time profile from tacho time signal
% [rpm,trpm] = tacho2rpm(x,fs,TLevel,Slope,PPR)
%
% rpm Instanteneous rpm
% t Time axis for rpm
%
% x Tacho time signal
% fs Sampling frequency (in Hz) for x
% TLevel Trigger level
% Slope +1, or -1 for positive and negative slope, respectively
% PPR Number of pulses per revolution for tacho signal
% NewFs (optional) sampling frequency for rpm. Default is fs
% FilterL Length of smoothing filter, Default is 5 pulses
%
% See also
%
% Copyright (c) 2009-2011 by Anders Brandt
% Email: abra@iti.sdu.dk
% Version: 1.0 2011-06-23
% 1.1 2012-01-15 Modified to work in Octave
% This file is part of ABRAVIBE Toolbox for NVA
% Parse inputs
if nargin == 5
NewFs = fs;
FilterL=5;
elseif nargin == 6
FilterL=5;
elseif nargin < 7
error('Wrong number of input parameters')
end
% Define time axis for tacho signal
t=makexaxis(x,1/fs);
% Get trigger times
xDiff=diff(sign(x-TLevel)); % Produces +/- 2 where trigger event
tDiff=t(2:end); % Diff shifts one sample
if Slope > 0
tTacho=tDiff(find(xDiff == 2)); % Tacho trigger positive slope instances
else
tTacho=tDiff(find(xDiff == -2)); % Tacho trigger negative slope instances
end
% Calculate rpm from time between tacho pulses. Assign rpm to second tacho
% pulse of each pair
rpmt=60/PPR./diff(tTacho); % Instantaneous rpm values
% Average each pair of rpm estimates and put on time of middle tacho pulse
rpmt=0.5*(rpmt(1:end-1)+rpmt(2:end));
tTacho=(tTacho(2:end-1)); % diff shifted one sample
% Smooth to obtain more stable values
if FilterL > 1
a=1;
b=1/FilterL*ones(1,FilterL);
if strcmp(checksw,'MATLAB')
rpmt=filtfilt(b,a,rpmt);
else % Octave filtfilt does not work
rpmt=filter(b,a,rpmt);
end
end
% Interpolate to NewFs resolution
t=(0:1/NewFs:t(end))';
rpm=interp1(tTacho,rpmt,t,'linear','extrap');
参考文献
[1] Shufeng Ai. Research on Order Tracking and Teager-Huang Transform Based Gear Crack Fault Diagnosis[C]. 2011 Second International Conference on Mechanic Automation and Control Engineering, Hohhot, 2011: 6041-6044.
[2] Z. Xiaofeng, S. Haibo and S. Wenli. An Implementation Method for the Transmission Fault Diagnosis Based on Order Tracking[C]. 2012 IEEE International Conference on Computer Science and Automation Engineering (CSAE), Zhangjiajie, 2012: 312-315.
[3] J. Hong, W. Guangrui, Z. Xiangfeng, et al. Application of Spectral Kurtosis and Vold-Kalman Filter Based Order Tracking in Wind Turbine Gearbox Fault Diagnosis[C]. 2017 9th International Conference on Modelling, Identification and Control (ICMIC), Kunming, 2017: 49-53.
