minphase.m

很多matlab的源代码

function [nn,dd]=minphase(ty,n,d)
% MINPHASE Convert Non-minmum Phase to Minimum Phase Transfer function
%
%	[NS,DS] = MINPHASE(TY,NU,DU) Finds a minimum phase transfer function 
%	TY = 's' for a CT HU(s) or TY = 'z' for a DT HU(z)
%	NU, DU are the num and den of the given transfer function.
%	NOTE: An allpass transfer function is used.
%	NS, DS are num and den of the minimum phase transfer function HS.
%
%	MINPHASE (with no input arguments) invokes the following example:
%          
%	% Convert H(z) =(z-4)/(z*z*z-5z*z+6*z) to a minimum phase TF 
%	  >>nu=[1 -4];du=[1 -5 6 0]
%	  >>rtn=roots(nu),rtd=roots(du)		%Check roots
%	  >>[ns,ds] = minphase('z',nu,du)	%Convert to min phase
%	  >>rtn=roots(ns),rtd=roots(ds)		%Check roots again


% 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 minphase,disp('Strike a key to see results of the example')
pause,nu=[1 -4];du=[1 -5 6 0];rtn=roots(nu),rtd=roots(du)
[nm,dm]=minphase('z',nu,du),rtn=roots(nm),rtd=roots(dm),return,end

while n(1)==0,n(1)=[];end,
while d(1)==0,d(1)=[];end,
g=n(1)/d(1);n=n/n(1);d=d/d(1);

%Min Phase Numerator
r=roots(n);
if ty=='z',i=find(abs(r)>1);else,i=find(real(r)>0);end
nn=n;if ~isempty(i)
ru=cplxpair(r(i));r(i)=[];
m=1;
for j=1:length(ru),
if ty=='z',r=[r;conj(1/ru(j))];m=m*ru(j);else,r=[r;-ru(j)];end
end
nn=real(poly(r));
mn=real(m);
if ty=='z',nn=nn*mn;end
end

%Min Phase Denominator
r=roots(d);
if ty=='z',i=find(abs(r)>1);else,i=find(real(r)>0);end
dd=d;if ~isempty(i)
ru=cplxpair(r(i));r(i)=[];
m=1;
for j=1:length(ru),
if ty=='z',r=[r;conj(1/ru(j))];m=m*ru(j);else,r=[r;-ru(j)];end
end
dd=real(poly(r)); 
md=real(m);if ty=='z',nn=nn/md;end
end

nn=g*nn;if nn(1)<0,nn=-nn;end,if dd(1)<0,dd=-dd;end