tfmin.m

很多matlab的源代码

function [pn,pd]=tfmin(nn,dd,ty,tol)
% TFMIN Minimum Phase Transfer Function from Magnitude squared function.
%
%	[N,D] = TFMIN(NUM,DEN,TY,TOL) Minimum phase TF from |H|^2.
%       If TY = 's', NUM, DEN = num and den of H(s=jw)^2 = NUM(w)/DEN(w)
%	If TY = 'z', NUM, DEN = num and den of H(z) with SYMMETRIC coefficients
%       TY = 's', for H(s), TY = 'z' for H(z) [Default: TY = 's'].
%	TOL is the tolerance for repeated roots. [default: TOL = 0.001]
%       N, D are the numerator and denominator of min-phase TF
%
%	NOTE: The DC gains of Hmin and Horiginal are matched (only at dc) 
%
%       USAGE: For H(z)  = (2z*z+5z+2)/(4z*z+17z+4) with symmetric coeffs, use
%		>>ytm = tfmin([2,5,2],[4,17,4],'z')
%
%       TFMIN (with no input arguments) invokes the following example:
%
%       % Find the Min phase TF of a 3rd order BW filter for which
%	% |H(w)|^2 = 1/(1+w^6)
%         >>[nmin,dmin] = tfmin(1,[1 0 0 0 0 0 1],'s')


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998



if nargin==0,help tfmin,disp('Strike a key to see results of first example')
pause,[nmin,dmin]=tfmin(1,[1 0 0 0 0 0 1],'s'),return,end

if nargin<4,tol=0.001;end,if nargin<3,ty='s';end,ty=ty(1);
while nn(1)==0,nn(1)=[];end,while dd(1)==0,dd(1)=[];end
n=length(nn);d=length(dd);
if ty=='s'
if rem(n,2)==0,error('Num must have odd length'),return,end
if rem(d,2)==0,error('Den must have odd length'),return,end
if n>1,l=nn(2:2:n);if any(l),error('Num not function of w^2'),return,end,end
l=nn(1:2:n);if l(1)<0,l=-l;end
if any(any(l<0)),error('All Num coeffs must have same sign'),return,end
if d>1,m=dd(2:2:d);if any(m),error('Den not function of w^2'),return,end,end
m=dd(1:2:d);if m(1)<0,m=-m;end
if any(any(m<0)),error('All Den coeffs must have same sign'),return,end
elseif ty=='z'
f=nn-fliplr(nn);if any(f),error('Num not symmetric'),return,end
f=dd-fliplr(dd);if any(f),error('Den not symmetric'),return,end
else,error('Unrecognized type: use s or z'),return
end

c=0;if ty=='s',if dd(d)~=0,c=nn(n)/dd(d);c=sqrt(abs(c));end,end
if ty=='z',ds=sum(dd);if ds~=0,c=sum(nn)/ds;c=sqrt(abs(c));end,end

if ty=='s',nn=nn.*((-1).^((n-1:-1:0)/2));dd=dd.*((-1).^((d-1:-1:0)/2));end
rn=roots(nn);rd=roots(dd);
j=sqrt(-1);qr=real(rn);qi=imag(rn);
r=abs(qr-round(qr))<tol;qr=round(qr).*r+qr.*(1-r);
r=abs(qi-round(qi))<tol;qi=round(qi).*r+qi.*(1-r);
rn=cplxpair(qr+j*qi);
qr=real(rd);qi=imag(rd);
r=abs(qr-round(qr))<tol;qr=round(qr).*r+qr.*(1-r);
r=abs(qi-round(qi))<tol;qi=round(qi).*r+qi.*(1-r);
rd=cplxpair(qr+j*qi);

if ty=='z'
if ~isempty(rn),i=find(imag(rn)~=0);
if ~isempty(i),pi=rn(i);rn(i)=[];else,pi=[];end,pr=rn;
r1=rem(length(pr),2);r2=rem(length(pi),4);rtn=[r1;r2];
if any(rtn),error('Num is not reciprocal cojnugate symmetric'),return,end,end
if ~isempty(rd),i=find(imag(rd)~=0);
if ~isempty(i),qi=rd(i);rd(i)=[];else,qi=[];end,qr=rd;
r3=rem(length(qr),2);r4=rem(length(qi),4);rtd=[r3;r4];
if any(rtd),error('Den is not reciprocal cojnugate symmetric'),return,end,end
end
if ty=='s'
x=[];i=find(rn==0);if ~isempty(i);r1=rn(i);l=length(i)/2;
x=[x;rn(1:l)];rn(i)=[];end
y=[];i=find(rd==0);if ~isempty(i);r1=rd(i);l=length(i)/2;
y=[y;rd(1:l)];rd(i)=[];end
i=find(abs(real(rn))<tol);pp=cplxpair(rn(i));if ~isempty(i),rn(i)=[];end
while ~isempty(pp),if abs(pp(1)-pp(3))<tol,x=[x;pp(1:2)];pp(1:4)=[];end,end
i=find(abs(real(rd))<tol);pp=cplxpair(rd(i));if ~isempty(i),rd(i)=[];end
while ~isempty(pp),if abs(pp(1)-pp(3))<tol,y=[y;pp(1:2)];pp(1:4)=[];end,end
i=find(real(rn)<0);x=[x;rn(i)];i=find(real(rd)<0);y=[y;rd(i)]; %Find LHP roots
else
x=[];i=find(pr==1);if ~isempty(i);l=length(i)/2;x=[x;pr(1:l)];pr(i)=[];end
i=find(pr==-1);if ~isempty(i);l=length(i)/2;x=[x;pr(1:l)];pr(i)=[];end
i=find(abs(pr)<1);x=[x;pr(i)];
i=find(abs(abs(pi)-1)<tol);pp=cplxpair(pi(i));if ~isempty(i),pi(i)=[];end
while ~isempty(pp),if abs(pp(1)-pp(3))<tol,x=[x;pp(1:2)];pp(1:4)=[];end,end
i=find(abs(pi)<1);x=[x;pi(i)];
y=[];i=find(qr==1);if ~isempty(i);l=length(i)/2;x=[x;qr(1:l)];qr(i)=[];end
i=find(qr==-1);if ~isempty(i);l=length(i)/2;x=[x;qr(1:l)];qr(i)=[];end
i=find(abs(qr)<1);y=[y;qr(i)];
i=find(abs(abs(qi)-1)<tol);pp=cplxpair(qi(i));if ~isempty(i),qi(i)=[];end
while ~isempty(pp),if abs(abs(pp(1))-abs(pp(3)))<=tol,y=[y;pp(1:2)];pp(1:4)=[];
end,end
i=find(abs(qi)<1);y=[y;qi(i)];
end
pn=real(poly(x));pd=real(poly(y));
%lx=length(pn);ly=length(pd);ld=max(lx,ly);
%r=abs(pn-round(pn))<20*eps;pn=round(pn).*r+pn.*(1-r);
%r=abs(pd-round(pd))<20*eps;pd=round(pd).*r+pd.*(1-r);
%pn=[zeros(1,ld-lx) pn];pd=[zeros(1,ld-ly) pd];
%if c~=0,if ty=='s',pn=pn*c*pd(ld)/pn(ld);else,pn=pn*c*sum(pd)/sum(pn);end,end
%y=[pn;pd];
ln=length(pn);ld=length(pd);
if c~=0,if ty=='s',pn=pn*c*pd(ld)/pn(ln);else,pn=pn*c*sum(pd)/sum(pn);end,end