spec_fac.m

来自「matlab的大量实例代码」· M 代码 · 共 22 行

function [B,A] =spec_fac(num,den)
m=length(den)-1;
NN=conv(num,(-1).^[length(num)-1:-1:0].*num);
X=NN(1:2:end)';  X=0.5*[zeros(m-length(X),1); X];
p=roots(den); ii=find(p>0); p(ii)=-p(ii); A=poly(p); 
k=0; 
if m>1
   Xx = [(-1)^(m-1)*A(2) (-1)^m*A(1) zeros(1,m-2)];
else
   Xx=[A(1)];
end
V0 = Xx;
for i=2:m
   V0 = [0 0 V0(1:m-2)]; k = k+2;
   if k < m+1,
      V0(2) = (-1)^m*A(k+1);
      if k < m, V0(1) = (-1)^(m-1)*A(k+2); end
   end
   Xx = [Xx; V0];
end
B=[inv(Xx)*X]';