www.gusucode.com > signal 工具箱matlab源码程序 > signal/arburg.m
function [Aout,Eout,Kout] = arburg(x, p) %ARBURG AR parameter estimation via Burg method. % A = ARBURG(X,ORDER) returns the coefficients of the autoregressive (AR) % parametric signal model estimate of X using Burg's method. The model % has order ORDER, and the output array A has ORDER+1 columns. The % coefficients along the Nth row of A model the Nth column of X. If X is % a vector then A is a row vector. % % [A,E] = ARBURG(...) returns the final prediction error E (the variance % estimate of the white noise input to the AR model). % % [A,E,K] = ARBURG(...) returns the reflection coefficients (parcor % coefficients) in each column of K. % % % Example: % % Estimate input noise variance for AR(4) model. % A=[1 -2.7607 3.8106 -2.6535 0.9238]; % % % Generate noise standard deviations % % Seed random number generator for reproducible results % rng default; % noise_stdz=rand(1,50)+0.5; % % % Generate column vectors that have corresponding standard deviation % x = bsxfun(@times,noise_stdz,randn(1024,50)); % % % filter each column using the AR model. % y = filter(1,A,x); % % % Compute the estimated coefficients and deviations for each column % [ar_coeffs,NoiseVariance]=arburg(y,4); % % %Display the mean value of each estimated polynomial coefficient % estimatedA = mean(ar_coeffs) % % %Compare actual vs. estimated variances % plot(noise_stdz.^2,NoiseVariance,'k*'); % xlabel('Input Noise Variance'); % ylabel('Estimated Noise Variance'); % % See also PBURG, ARMCOV, ARCOV, ARYULE, LPC, PRONY, FILLGAPS. % Ref: S. Kay, MODERN SPECTRAL ESTIMATION, % Prentice-Hall, 1988, Chapter 7 % S. Orfanidis, OPTIMUM SIGNAL PROCESSING, 2nd Ed. % Macmillan, 1988, Chapter 5 % Copyright 1988-2014 The MathWorks, Inc. narginchk(2,2) % Checks if X is valid data signal.internal.sigcheckfloattype(x,'','arburg','X'); % convert to column vector before 2d check is performed if isvector(x) x = x(:); end validateattributes(x,{'numeric'},{'nonempty','finite','2d'},'arburg','X'); validateattributes(p,{'numeric'},{'positive','integer','scalar'},'arburg','ORDER'); % Cast to enforce precision rules p = double(p); if issparse(x), error(message('signal:arburg:Sparse')) end if size(x,1) < p+1 if isvector(x) error(message('signal:arburg:VectorTooSmall', p + 1)); else error(message('signal:arburg:MatrixTooSmall', p + 1)); end end N = size(x,1); Nchan = size(x,2); % By convention all polynomials are row vectors Aout = zeros(Nchan, p+1, class(x)); % Optional outputs Eout = zeros(1, Nchan, class(x)); Kout = zeros(p, Nchan, class(x)); for ichan = 1:Nchan % Initialization efp = x(2:end,ichan); ebp = x(1:end-1,ichan); % Initial error E = x(:,ichan)'*x(:,ichan)./N; a = zeros(1,p+1,class(x)); a(1) = 1; for m=1:p % Calculate the next order reflection (parcor) coefficient k = (-2.*ebp'*efp) ./ (efp'*efp + ebp'*ebp); Kout(m,ichan) = k; % Update the forward and backward prediction errors ef = efp(2:end) + k .* ebp(2:end); ebp = ebp(1:end-1) + k' .* efp(1:end-1); efp = ef; % Update the AR coeff. a(2:m+1) = a(2:m+1) + k .* conj(a(m:-1:1)); % Update the prediction error E = (1 - k'.*k)*E; end % Assign outputs Aout(ichan,:) = a; Eout(ichan) = E; end