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

    function varargout = pmem( xR, p, nfft, Fs, flag )
%PMEM   Power Spectrum estimate via MEM (Maximum Entropy Method).
%   PMEM has been replaced by PYULEAR.  PMEM still works but may be
%   removed in the future.  Type help PYULEAR for details.
%
%   See also PBURG, PCOV, PMCOV, PMTM, PMUSIC, PEIG, PWELCH, PERIODOGRAM.

%   Ref: S.L. Marple, DIGITAL SPECTRAL ANALYSIS WITH APPLICATIONS,
%              Prentice-Hall, 1987, Chapter 7

%   Author(s): J. McClellan, 9-17-95
%   Copyright 1988-2012 The MathWorks, Inc.
%     

narginchk(2,5);

[Mx,Nx] = size(xR);
xIsMatrix = Mx>1 & Nx>1;

if nargin < 1
   error(message('signal:pmem:Nargchk'))
elseif  isempty(p)
   error(message('signal:pmem:Empty'))
end
if issparse(xR)
   error(message('signal:pmem:Sparse'))
end
if nargin < 5,   flag = [];   end
if nargin < 4,   Fs = [];     end  
if nargin < 3,   nfft = [];   end
if ischar(nfft)
   tflag = nfft; nfft = Fs; Fs = flag; flag = tflag;
elseif ischar(Fs)
   tflag = Fs;  Fs = flag; flag = tflag;
end

if isempty(nfft),    nfft = 256;  end
if isempty(Fs),      Fs = 2;      end
corr_flag = 0;     %<---- not expecting correlation
if ~isempty(flag)
   flag = upper(flag);
   if  (~isempty( strfind(flag,'CORR') )),  corr_flag = 1;  end
end

if( corr_flag )   %-- might be correlation matrix
   if Mx~=Nx
      error(message('signal:pmem:MatrixMustBeSquare'))
   elseif  norm(xR'-xR) > 100*eps
      error(message('signal:pmem:SignalErr'))
   end
end

if  ~xIsMatrix
   [RR,lags] = xcorr(xR,p,'biased'); %#ok
   xR = toeplitz( RR(p+1:2*p+1), RR(p+1:-1:1) );
else  %-- Matrix case
   if  p>= Mx
      error(message('signal:pmem:ColLengthMustBeGtOrder'))
   elseif  ~corr_flag
      xR = corrcoef( xR.' );
   end
end
a = [1;  xR(2:p+1,2:p+1)\(-xR(2:p+1,1)) ];
EE = abs( xR(1,1:p+1) * a );

Spec2 = abs( fft( a, nfft ) ) .^ 2;
Pxx = EE*ones(size(Spec2)) ./ Spec2;

%--- Select first half only, when input is real
if isreal(xR)   
    if rem(nfft,2),    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
else
    select = (1:nfft)';
end

Pxx = Pxx(select);
ff = (select - 1)*Fs/nfft;

if nargout == 0
   newplot;
   plot(ff,10*log10(Pxx)), grid on
   xlabel('Frequency'), ylabel('Power Spectrum Magnitude (dB)');
   title('Maximum Entropy Spectral Estimate')
end

if nargout >= 1
    varargout{1} = Pxx;
end
if nargout >= 2
    varargout{2} = ff;
end
if nargout >= 3
    varargout{3} = a;
end