freq2levy.m

来自「系统辨识的M文件」· M 代码 · 共 36 行

function [num,den]=FREQ2LEVY(H,w,na,nb)
% FREQ2LEVY Analog filter weighted least squares fit to frequency response data 
%   by Levy algorithm..
%   [NUM,DEN]=FREQ2LEVY(H,W,NA,NB) gives real numerator and denominator 
%   coefficients NUM and DEN of orders NB and NA respectively, where
%   H is the desired complex frequency response of the system at frequency
%   points W, and W contains the frequency values in radians/s.
%   FREQ2LEVY yields a filter with real coefficients.  This means that it is 
%   sufficient to specify positive frequencies only; the filter fits the data 
%   conj(H) at -W, ensuring the proper frequency domain symmetry for a real 
%   filter.

P=real(H);  Q=imag(H);  M=length(w);
for i=1:2:na+nb+1
    L(i)=sum(w.^(i-1));   S(i)=sum(w.^(i-1).*P);
    S(i+1)=sum(w.^i.*Q);  U(i+2)=sum(w.^(i+1).*(P.^2+Q.^2));
end
ii=[1,1]; 
for i=1:na
    ii=[ii  -ii(2*i-1:2*i)];
end
Gxx=L(1:na).*ii(1:na);
Gxy=S(2:nb+2).*ii(1:nb+1);
Gyy=U(3:nb+3).*ii(1:nb+1);
for i=1:na-1
    Gxx=[Gxx; -Gxx(1,2:na)    L(na+i)*ii(na-i)];
    Gxy=[Gxy; -Gxy(i,2:nb+1)  S(nb+2+i)*ii(na-i)];
end
for i=1:nb
    Gyy=[Gyy;  -Gyy(i,2:nb+1)  U(nb+i+3)*ii(nb+1-i)];
end
G=[Gxx,  Gxy; Gxy',   Gyy];
V=[S(1:na).*ii(1:na),  U(2:nb+2).*ii(1:nb+1)]';
V1=inv(G)*V;   num=V1(nb+1:-1:1)';   den=[V1(na+nb+1:-1:nb+2)' 1];
num=num/den(1);   den=den/den(1);