phdiff.m

PHDIFF函数是matlab函数

% PHDIFF  ( MatLinks) Finds the phase difference between two oscillations.
%
%    PHDIFF(X,Y,TS) finds the phase difference between oscillations X and Y,
%    representing two phase-shifted waveforms having the same frequency, for
%    the sampling time TS between successive data points in X and Y.  (X and Y
%    must be sampled at the same rate.)  This phase difference is equal to the
%    angular difference <X - <Y.  The default value of TS=1.
%
%    PHDIFF(...) by itself plots X, Y and their phase difference.
%
%    P=PHDIFF(...) returns the phase lag, in radians, of Y following X.  The
%    phase lag is thus always between 0 and 2*pi.  (In certain situations a
%    phase "lag" P>pi may be interpreted as a phase lead.)  P has two elements:
%    the first is the phase difference between the valleys of the two
%    oscillations, the second element is that of the peaks.  MEAN(P) thus gives
%    the mean phase difference.  It is useful to distinguish the two phase
%    difference for asymmetric periodic oscillations.  For symmetric oscilla-
%    tions, therefore, either element of P should be equivalent to the other
%    (P(1) = P(2)).
%
%    [P,W] = PHDIFF(...) also returns the mean frequency W, in rad/s, of the
%    two oscillations.  Thus P/W will convert the phase difference to the equi-
%    valent (asymmetric) time delay in seconds.
%
%    [P,W,R] = PHDIFF(...) also returns the mean amplitude ratio R = X/Y.  R
%    has two elements like P, the first being the amplitude ratio of the
%    oscillation valleys, the second being that of their peaks.
%
%    If X and Y have frequencies which are dissimilar by more than 2*pi*TS then
%    the return value for P=10 (chosen arbitrarily to be > 2*pi).
%
%    If the phase difference is undetectable then P=W=R = 0 is returned.
%
%    Note that small phase differences are estimated with relatively poor
%    precision when the number of oscillations or data point is small.
%
%    See also SETTLING, STEADY.
%
%    Type HELP MATLINKS for a full listing of all MatLinks ToolChest functions.
%
function [p, w, r] = phdiff(x, y, ts)
%===============================================================================
%  Copyright  1998,2000 Julian Andrew de Marchi, Ph.D. ([e-mail])
%  Use & distribution covered by GNU General Public License (www.gnu.org)
%===============================================================================

%------------------
% parse the inputs
%------------------
if (nargin == 0),  error('No oscillation vector X supplied.');
elseif (nargin == 1),  error('No oscillation vector Y supplied.');
elseif (nargin < 3),  ts = 1;
end;

x = x(:);  y = y(:);
if (length(x) ~= length(y)),  error('Oscillations X and Y must have the same length.');
elseif (ts <= 0),  error('Time interval TS must be a positive number.');  end;

%--------------------------------------------------------------------------
% find the relative peak locations
%--------------------------------------------------------------------------
xs = steady(x);  ys = steady(y);
if (xs & ys),
  ss = max([xs ys]);  xx = x(ss:length(x));  yy = y(ss:length(y));
else
  xs = settling(x);  ys = settling(y);
  if (xs & ys),
    ss = min([xs ys]);  xx = x(1:ss);  yy = y(1:ss);  ss = 1;
  else
    warning('Warning: X and Y do not settle enough to estimate their phase difference.')
    p = [0 0];  w = 0;  r = [0 0];  return;
  end;
end;

%[xv xp] = findpeak(xx,max(abs(xx)/10));  [yv yp] = findpeak(yy,max(abs(yy)/10));
[xv xp] = findpeak(xx);
if (xv(1) == 1)
  xv = xv(2:length(xv));
elseif (xp(1) == 1)
  xp = xp(2:length(xp));
end;

[yv yp] = findpeak(yy);
if (yv(1) == 1)
  yv = yv(2:length(yv));
elseif (xp(1) == 1)
  yp = yp(2:length(yp));
end;

if (xv(1) == 0 | yv(1) == 0),
  p = [0 0];  w = 0;  r = [0 0];  return;
elseif (length(xv) > length(yv)),
  xv = xv(1:length(yv));
elseif (length(yv) > length(xv)),
  yv = yv(1:length(xv));
end;
if (xp(1) == 0 | yp(1) == 0),
  p = [0 0];  w = 0;  r = [0 0];  return;
elseif (length(xp) > length(yp)),
  xp = xp(1:length(yp));
elseif (length(yp) > length(xp)),
  yp = yp(1:length(xp));
end;
d0v = mean(yv - xv);  d0p = mean(yp - xp);  d0 = [d0v d0p];
%r = mean((x(xvp) - (min(x(xvp))+max(x(xvp)))/2) ./ (y(yvp) - (min(y(yvp))+max(y(yvp)))/2));
rv = mean((x(xv) - (min(x(xv))+max(x(xp)))/2) ./ (y(yv) - (min(y(yv))+max(y(yp)))/2));
rp = mean((x(xp) - (min(x(xv))+max(x(xp)))/2) ./ (y(yp) - (min(y(yv))+max(y(yp)))/2));
r = [rv rp];

%-------------------------------------------------
% compute the mean frequency and phase difference
%-------------------------------------------------
w = 2*pi / (ts*abs(mean([diff(xv) diff(xp) diff(yv) diff(yp)])));
p = d0*w*ts;
p(find(p<0)) = p(find(p<0)) + 2*pi;

%-------------------------------------------------------
% plot the peak if there's no output variable
%-------------------------------------------------------
if (nargout==0),
  hold off, plot(ts*(1:length(x)), [x y]), hold on,
  plot(ts*(ss+[xv xp]-1), xx([xv xp]), 'go'), plot(ts*(ss+[yv yp]-1), yy([yv yp]), 'ro');
  if (nargin==3), xlabel('time (s)');
  else xlabel('time (samples)'); end;
  ylabel('X and Y'), title(['Phase diff.=' num2str(mean(p)) '+/-' num2str(abs(diff(p))) ' rad. (' num2str(mean(p)*180/pi) '+/-' num2str(abs(diff(p))*180/pi) ' deg.) at ' num2str(w) ' rad/s (' num2str(w/(2*pi)) ' Hz.)']), zoom on;
  p = [p; w 0; r];
end;


%===============================================================================
% End-of-File
%===============================================================================