dir21adr.m

来自「关于数字信号处理的一些matlab例程」· M 代码 · 共 17 行

function [K,C]=dir21adr(b,a)
a1=a(1);a=a/a1;b=b/a1;
M=length(b);N=length(a);
if M>N
   error('  *** length of b must be <=length of a ***')
end
b=[b,zeros(1,N-M)];K=zeros(1,N-1);
A=zeros(N-1,N-1);C=b;
for m=N-1:-1:1
A(m,1:m)=-a(2:m+1)*C(m+1);
K(m)=a(m+1);
J=fliplr(a);
a=(a-K(m)*J)/(1-K(m)*K(m));
a=a(1:m);
C(m)=b(m)+sum(diag(A(m:N-1,1:N-m)));
end