www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/trench.m
function [amat,cmat,pf,gamf,gamb] = trench(c,r) %TRENCH Trench recursion (non-symmetric Toeplitz matrix) % [amat,cmat, pfe,gamf,gamb] = trench(c,r) % c - first column of Toeplitz matrix, r(k), k=0,...,M % r - first row of Toeplitz matrix, r(-k), k=0,...,M % defaults to c'. % amat - estimated forward AR predictors, of orders 1 through M % the k-th column contains the AR(k) coefficients % cmat - estimated backward AR predictors, of orders 1 through M % the k-th column contains the AR(k) coefficients % pfe - final prediction error variance for orders 1 through M % the k-th element contains the value for order k % gamf - estimated forward reflection coefficients % gamb - estimated backward reflection coefficients % % Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved. % $Revision: 1.5 $ % A. Swami January 20, 1995 % RESTRICTED RIGHTS LEGEND % Use, duplication, or disclosure by the Government is subject to % restrictions as set forth in subparagraph (c) (1) (ii) of the % Rights in Technical Data and Computer Software clause of DFARS % 252.227-7013. % Manufacturer: United Signals & Systems, Inc., P.O. Box 2374, % Culver City, California 90231. % % This material may be reproduced by or for the U.S. Government pursuant % to the copyright license under the clause at DFARS 252.227-7013. % --------------------------------------------------------- c = c(:); if (exist('r') ~= 1) r = c'; end if (c(1) ~= r(1)) disp(['r(1) ~= c(1); setting r(1) = c(1) ']); r(1) = c(1); end m = min(length(c), length(r)) - 1; c = c(1:m+1); r = r(1:m+1); r = r(:); pf = zeros(m+1,1); df = zeros(m,1); db = zeros(m,1); gamf = zeros(m,1); gamb = zeros(m,1); pf(1) = c(1); df(1) = c(2); db(1) = r(2); avec = [1]; cvec = [1]; amat = zeros(m+1,m); cmat = zeros(m+1,m); for k=1:m gamf(k) = -df(k) / pf(k); gamb(k) = -db(k) / pf(k); if ( gamf(k)*gamb(k) > 1) disp(['gamf * gamb > 1 at order',int2str(k)]) disp(['stopping iterations']) break end pf(k+1) = pf(k) * (1 - gamf(k)*gamb(k)); aold = avec; avec = [avec; 0] + gamf(k) * [0; cvec]; cvec = [0; cvec] + gamb(k) * [aold ;0]; if (k < m) df(k+1) = c(k+2:-1:2).' * avec; db(k+1) = r(2:1:k+2).' * cvec; end amat(1:k+1,k) = avec; cmat(m+1-k:m+1,k) = cvec; end pf = pf(2:m+1); return