gusucode.com > 信号处理工具箱 - signal源码程序 > signal\signal\signal\rlevinson.m
function [R,U,kr,e]=rlevinson(a,efinal) %RLEVINSON Reverse Levinson-Durbin Recursion. % R=RLEVINSON(A,Efinal) computes the autocorrelation coefficients, R based % on the prediction polynomial A and the final prediction error Efinal, % using the stepdown algorithm. A should be a minimum phase polynomial % and A(1) is assumed to be unity. % % [R,U]=RLEVINSON(...) returns a (P+1)x(P+1) upper triangular matrix, U, % that holds the i'th order prediction polynomials Ai, i=1:P, where P % is the order of the input polynomial, A. % % The format of this matrix U, is: % % [ 1 a1(1)* a2(2)* ..... aP(P) * ] % [ 0 1 a2(1)* ..... aP(P-1)* ] % U = [ .................................] % [ 0 0 0 ..... 1 ] % % from which the i'th order prediction polynomial can be extracted % using Ai=U(i+1:-1:1,i+1)'. The first row of U contains the % conjugates of the reflection coefficients, and the K's may be % extracted using, K=conj(U(1,2:end)). % % [R,U,K]=RLEVINSON(...) returns the reflection coefficients in K. % [R,U,K,E]=RLEVINSON(...) returns the prediction error of descending % orders P,P-1,...,1 in the vector E. % % See also LEVINSON. % References: [1] S. Kay, Modern Spectral Estimation, % Prentice Hall, N.J., 1987, Chapter 6. % [2] P. Stoica R. Moses, Introduction to Spectral Analysis % Prentice Hall, N.J., 1997, Chapter 3. % % Author(s): A. Ramasubramanian % Copyright (c) 1988-98 by The MathWorks, Inc. % $Revision: 1.6 $ $Date: 1998/08/27 13:17:50 $ % Some preliminaries first a = a(:); % Convert to a column vector if not already so if a(1)~=1, warning(['First coefficient of the prediction polynomial ' ... 'was not unity.']); a = a/a(1); end p = length(a); if p < 2 error('Polynomial should have at least two coefficients'); end U = zeros(p,p); % This matrix will have the prediction polynomials % of orders 1:p U(:,p) = conj(a(end:-1:1)); % Prediction coefficients of order p p = p-1; % First we find the prediction coefficients of smaller orders and form the % Matrix U % Initialize the step down e(p) = efinal; % Prediction error of order p % Step down for k=p-1:-1:1 [a,e(k)] = levdown(a,e(k+1)); U(:,k+1) = [conj(a(end:-1:1).');zeros(p-k,1)]; end e0 = e(1)/(1-abs(a(2)^2)); % Because a[1]=1 (true polynomial) U(1,1) = 1; % Prediction coefficient of zeroth order kr = conj(U(1,2:end)); % The reflection coefficients kr = kr.'; % To make it into a column vector % Once we have the matrix U and the prediction error at various orders, we can % use this information to find the autocorrelation coefficients. % Initialize recursion k = 1; R0 = e0; % To take care of the zero indexing problem R(1) = -conj(U(1,2))*R0; % R[1]=-a1[1]*R[0] % Actual recursion for k = 2:1:p R(k) = -sum(conj(U(k-1:-1:1,k)).*R(end:-1:1).')-kr(k)*e(k-1); end % Include R(0) and make it a column vector. Note the dot transpose R = [R0 R].'; % [EOF] rlevinson.m