www.gusucode.com > signal 工具箱matlab源码程序 > signal/tachorpm.m

    function varargout = tachorpm(x,fs,varargin)
%TACHORPM Extract an RPM signal from tachometer pulses
%   RPM = TACHORPM(X,Fs) extracts a rotational speed signal vector, RPM,
%   from a tachometer pulse signal vector, X, with sampling frequency, Fs,
%   in hertz. RPM has the same number of elements as X.
%
%   [RPM,T] = TACHORPM(X,Fs) returns the RPM signal time vector, T,
%   measured in seconds.
%
%   [RPM,T,TP] = TACHORPM(X,Fs) returns a vector of detected pulse
%   times, TP, measured in seconds.
%
%   [...] = TACHORPM(...,'PulsesPerRev',PPR) specifies the number of
%   tachometer pulses per revolution, PPR. If PPR is not specified, it
%   defaults to 1.
%
%   [...] = TACHORPM(...,'StateLevels',SL) sets the levels used to identify
%   pulses. Specify SL as a two-element real vector with first and second
%   elements corresponding to the lower and upper state levels of the input
%   waveform. Choose the state levels so that all pulse edges cross within
%   10% of both of them. If SL is not specified, the levels are computed
%   automatically using the histogram of the input waveform.
%
%   [...] = TACHORPM(...,'OutputFs',OFs) specifies the sampling frequency
%   OFs, of the outputs RPM and T. If OFS is not specified, it defaults to
%   the tachometer sampling frequency, Fs.
%
%   [...] = TACHORPM(...,'FitType',FT) uses the fitting method
%   specified by FT to generate RPM. FT can be
%   'smooth' or 'linear':
%     'smooth' - fit a least-squares B-spline to pulse rpm values
%     'linear' - interpolate linearly between pulse rpm values
%   If FT is not specified, it defaults to 'smooth'.
%
%   [...] = TACHORPM(...,'FitPoints',FP) uses FP breakpoints in the
%   least-squares B-spline. The number of points, FP, is a trade-off
%   between curve smoothness (fewer points) and closeness to the underlying
%   data (more points). Choosing too many points can result in overfitting
%   the data. If unspecified, FP defaults to 10. If FT is
%   'linear', FP is ignored.
%
%   TACHORPM(...) with no output arguments creates plots containing the
%   generated rpm signal, the tachometer signal, and the detected pulses.
%
%   % EXAMPLE:
%   %   Load a tachometer signal.
%   load tacho.mat
%
%   % Compute and visualize the rpm signal using a large number of fit
%   % points to capture the rpm peak.
% 
%   tachorpm(Yn,fs,'FitPoints',100)
%   xlim([0 .1])
%   ylim([500 900])
%
%   See also RPMORDERMAP, ORDERTRACK, ORDERWAVEFORM, STATELEVELS

%   References:
%     [1] Brandt, Anders. Noise and Vibration Analysis: Signal Analysis 
%         and Experimental Procedures. Chichester, UK: John Wiley & Sons,
%         2011.

% Copyright 2015-2016 MathWorks, Inc.

narginchk(1,8);
nargoutchk(0,3);

% Parse inputs
[ppr,ofs,ft,sl,fp] = parseInputs(fs,varargin{:});

validateInputs(x,fs,ppr,ofs,sl,fp);

% Cast to enforce precision rules
% Since risetime, falltime, and statelevels won't accept single, cast to
% double and recast the output to x_class.
x_class = class(x);
x = double(x);
fs = double(fs);
ppr = double(ppr);
sl = double(sl);
ofs = double(ofs);

x = x(:);

% Use statelevels if levels have not been provided
if isempty(sl)
  sl = statelevels(x);
  sl = sl(:);
end

% Construct time vector
t = (0:length(x)-1)'/fs;

% Compute the rising and falling edges of the tachometer pulses using
% risetime and falltime.
[~,LT1] = risetime(x,fs,'StateLevels',sl','Tolerance',10,...
  'PercentReferenceLevels',[50 51]);
[~,LT2] = falltime(x,fs,'StateLevels',sl','Tolerance',10,...
'PercentReferenceLevels',[50 51]);

% Error out if no edges were detected
if isempty(LT1) || isempty(LT2)
  error(message('signal:rpmmap:StateLevelsEmpty'));
end

% Make sure number of rising and falling edges is the same. If not, remove
% extra from the end of the signal
if length(LT1)~=length(LT2) 
  l = min([length(LT1) length(LT2)]);
  LT1 = LT1(1:l);
  LT2 = LT2(1:l);
end

% Compute the times of the center of each pulse
tCenter= mean([LT1(:) LT2(:)],2);

% Determine the period between pulse centers
T = diff(tCenter); 

% Calculate RPM values at the center points between pulse centers. 
rpmPulse = 60/ppr./T; 
tPulse = interp1(1:length(tCenter),tCenter,1.5:length(tCenter)-0.5);

% Compute the output rpm signal. Use either least-squared splines or linear
% interpolation.
tout = interp1(0:length(t)-1,t,0:fs/ofs:length(t)-fs/ofs,'linear','extrap')';

if strcmp(ft,'smooth')
  % Use a least-squares spline to produce the rpm signal
  b = linspace(tPulse(1),tPulse(end),fp);
  sp = spline(b,rpmPulse(:)'/spline(b,eye(length(b)),tPulse(:)'));
  rpm = ppval(sp,tout);
else
  % Use linear interpolation to produce the rpm signal
  rpm = interp1(tPulse,rpmPulse,tout,'linear','extrap');
end

% Cast output if x was single
if isequal(x_class,'single')
  rpm = single(rpm);
  tout = single(tout);
  tCenter = single(tCenter);
end

if nargout == 0 
   plotTacho(t,tout,tPulse,tCenter,x,rpm,rpmPulse,sl);
elseif nargout == 1
  varargout{1} = rpm(:);
elseif nargout == 2
  varargout{1} = rpm(:);
  varargout{2} = tout(:);
elseif nargout == 3
  varargout{1} = rpm(:);
  varargout{2} = tout(:);
  varargout{3} = tCenter(:);
end
%--------------------------------------------------------------------------
function [ppr,ofs,ft,sl,fp] = parseInputs(fs,varargin)

% Default values
ppr = 1;
ofs = fs;
ft = 'smooth';
sl = [];
fp = 10;

% Check that name-value inputs come in pairs and are all strings
if rem(numel(varargin),2)
  error(message('signal:rpmmap:NVMustBeEven'));
end

% Parse Name-value pairs
p = inputParser;
p.addParameter('PulsesPerRev',ppr);
p.addParameter('OutputFs',ofs);
p.addParameter('FitType',ft);
p.addParameter('StateLevels',sl);
p.addParameter('FitPoints',fp);
parse(p,varargin{:});
ppr = p.Results.PulsesPerRev;
ofs = p.Results.OutputFs;
ft = p.Results.FitType;
sl = p.Results.StateLevels;
fp = p.Results.FitPoints;

% Check that ft is a match to allowed values
ft = validatestring(ft,{'linear','smooth'},'tachorpm','FT'); 

% Make sl a column vector
if ~isempty(sl)
  sl = sl(:);
end
%--------------------------------------------------------------------------
function validateInputs(x,fs,ppr,ofs,sl,fp)
validateattributes(x,{'numeric'},...
  {'real','finite','nonnan','vector'},'tachorpm','X');
validateattributes(fs,{'numeric'},...
  {'real','positive','nonsparse','finite','nonnan','scalar'},'tachorpm','Fs');
validateattributes(ppr,{'numeric'},...
  {'real','positive','nonsparse','finite','nonnan','scalar'},'tachorpm','PPR');
validateattributes(ofs,{'numeric'},...
  {'real','positive','nonsparse','finite','nonnan','scalar'},'tachorpm','OFs');
validateattributes(fp,{'numeric'},...
  {'real','positive','nonsparse','finite','nonnan','scalar','>=',2},'tachorpm','FP');
if ~isempty(sl)
  validateattributes(sl,{'numeric'},...
    {'real','finite','nonsparse','nonnan','vector','numel',2},'tachorpm','SL');
end
%--------------------------------------------------------------------------
function plotTacho(t,tout,tPulse,tCenter,x,rpm,rpmPulse,sl)
  [~,E,U]=engunits(t,'unicode','time');
  xlab = [getString(message('signal:spectrogram:Time')) ' (' U ')'];
  
  % Create subplots, change their size, and make them invisible until the
  % plot command is called (which makes them visible automatically).
  p1 = subplot(2,1,1);
  p2 = subplot(2,1,2);
  linkaxes([p1,p2],'x');

  % Make the rpm plot larger and the pulse plot larger
  p1p = get(p1,'position');
  p2p = get(p2,'position');
  p1p(2) = 0.3612;
  p1p(4) = 0.55;
  p2p(4) = 0.175;
  set(p1,'position',p1p)
  set(p2,'position',p2p)

  % RPM plot
  axes(p1);
  hold on
  plot(p1,tout*E,rpm); 
  plot(p1,tPulse*E,rpmPulse,'+');
  ylabel(getString(message('signal:rpmmapplot:RPMs')));
  legend(getString(message('signal:rpmmapplot:rpmsignal')),...
    getString(message('signal:rpmmapplot:rpmpulses')),'Location','SouthEast');
  grid on;
  title(getString(message('signal:rpmmapplot:rpmSignal')));
  set(p1,'xTickLabels',[])
  axis('tight');
  hold off;

  % Pulse and tacho plot
  axes(p2);
  hold on
  plot(p2,tCenter*E,mean(interp1(t,x,tCenter))*ones(size(tCenter)),'+', ...
    'Color',[0.8500 0.3250 0.0980]);
  plot(p2,tout*E,(sl(1)) * ones(size(tout)),'--','Color', ...
    [0.4660 0.6740 0.1880],'MarkerSize',1);
  plot(p2,tout*E,(sl(2)) * ones(size(tout)),'--','Color',...
    [0.4660 0.6740 0.1880],'MarkerSize',1);
  plot(p2,t*E,x,'Color',[0 0.4470 0.7410]);

  xlabel(xlab);
  ylabel(getString(message('signal:rpmmapplot:volts')));
  grid on;
  title(getString(message('signal:rpmmapplot:tachoSignal')));
  legend(getString(message('signal:rpmmapplot:detectedpulses')), ...
    getString(message('signal:rpmmapplot:statelevels')));
  axis('tight');
  hold off;

  % Make the plots tight in x and give a margin in y
  axes(p1);
  axis('tight');
  yl1 = get(p1,'ylim');
  yl2 = get(p2,'ylim');
  set(p2,'ylim',[-.1 .1]*abs(diff(yl2))+yl2);
  set(p1,'ylim',[-.1 .1]*abs(diff(yl1))+yl1);

  % Create tags
  p1.Tag = 'RPM';
  p2.Tag = 'Tacho';