beatfreq.m

BEATFREQ是matlab函数

% BEATFREQ  ( MatLinks) Finds the beat frequency of two oscillations.
%
%    BEATFREQ(X) finds the beat frequency of the data vector X.
%
%    BEATFREQ(X,W) ignores all multiple local (noise) peaks occurring within a
%    moving window of length W.  The default value for W=1.
%
%    BEATFREQ(...) by itself plots X and its beat with a dashed line.
%
%    Y=BEATFREQ(...) returns the beating functions within X.  Y has two columns,
%    each which contains one or the other of the constituent oscillating
%    functions causing the beat frequency phenomenon.
%
%    See also ENVELOPE.
%
%    Type HELP MATLINKS for a full listing of all MatLinks ToolChest functions.
%
function [y] = beatfreq(x, w)
%===============================================================================
%  Copyright  1998,2000 Julian Andrew de Marchi, Ph.D. ([e-mail])
%  Use & distribution covered by GNU General Public License (www.gnu.org)
% _____
% notes:
%   this code is buggy - JAD
%===============================================================================

%------------------
% parse the inputs
%------------------
if (nargin==0), error('No data vector X supplied.');
elseif (nargin<2), w=1;
elseif (w~=fix(w) | w<=0),
  error('The window length W must be a positive integer.');
end;

%-----------------------
% measure the envelopes
%-----------------------
[i,j] = peaks(x,0,w);
if (length(i)+length(j)<w+1),
  error('Not enough oscillation in X--envelope detection would be senseless.');
end;
len=min([length(i) length(j)]);  i=i(1:len);  j=j(1:len);
y = [spline(i,x(i),1:length(x))' spline(j,x(j),1:length(x))'];

%--------------------------
% cross the envelopes over
%--------------------------
[ii,jj]=peaks(abs(y(:,1)-y(:,2)));
inxi=1; inxj=1;
for inx=1:length(ii),
  while(i(inxi)<=ii(inx) & inxi<length(i)), inxi=inxi+1; end;
  while(j(inxj)<=ii(inx) & inxj<length(j)), inxj=inxj+1; end;
  if (inxi~=length(i) & inxj~=length(j)),
    dummyi=i(inxi:length(i)); dummyj=j(inxj:length(j));
    i=[i(1:inxi-1) dummyj]; j=[j(1:inxj-1) dummyi];
  end;
end;
y = [spline(i,x(i),1:length(x))' spline(j,x(j),1:length(x))'];

%---------------------------------------------
% plot the beat if there's no output variable
%---------------------------------------------
if (nargout==0),
  hold off,plot(x),hold on,plot(1:length(y),y(:,1),'g:',1:length(y),y(:,2),'c:'),
  xlabel('i'),ylabel('x(i)'),title(['Envelope (w=' num2str(w) ')']);
end;

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