lp2iir.m

很多matlab的源代码

function [nn,dd] = lp2iir(ty,ty1,n,d,sf,f,fc)
% LP2IIR Convert LPP to required IIR digital filter.
%
%       [NIIR,DIIR] = LP2IIR(TY,TY1,NLP,DLP,SF,F,FC) Convert LPP to IIR filter
%	TY = 'lp', 'hp', 'bp', or 'bs',  SF = sampling frequency in HERTZ 
%	TY1 = 'a' if conversion is from an analog LPP (the A2D transformation)  
%	TY1 = 'd' if conversion is from a digital LPP (the D2D transformation)
%	NLP, DLP are numerator/denominator of LPP in descending order
%	F contains the desired band edge(s) (Hz); (F = [F1 F2] for 'bp', 'bs')  
%	If TY1='d', FC = cutoff of the LPP (HZ) [Default: 1/(2*pi) Hz].
%	If TY1='a', FC=optional center frequency for BP and BS
%
%	NOTE: If TY1 = 'a', F must contain UNWARPED frequencies
%	NIIR, DIIR contain the coefficients of H(z) in descending order.
%
%	USAGE: Convert a 3rd order BW LPP to hp with cutoff = 100 Hz and 
% 	sf = 300Hz, using: >>[nh,dh] = lp2iir('hp','a',1,[1 2 2 1],300,100);
%
%	LP2IIR (with no input arguments) invokes the following example:
%
%        % Convert a DIGITAL LPF H(z)=z/(z-0.5) to BS with cutoffs [80 120]Hz,
%        % SF = 500Hz and assumed LPF cutoff=50Hz and plot the results
%        >>[nb,db] = lp2iir('bs','d',[1 0],[1 -0.5],500,[80 120],50); 
%        >>tfplot('z',nh,dh,[0 .5],0,1); 


% 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 lp2iir,disp('Strike a key to see results of last example')
pause,[nb,db]=lp2iir('bs','d',[1 0],[1 -0.5],500,[80 120],50),
tfplot('z',nb,db,[0 0.5],0,1);return,end

%%%%%%%%%%%%  USES new POLYMAP  %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%  ALSO USES DW (bandwidth+center freq) if given %%%%%%%


%if nargin<7,fc=1/2/pi;end
if length(f)==2,ws=sum(f)*pi/sf;wd=abs(diff(f))*pi/sf;end

if ty1=='d'
if nargin<7,fc=1/2/pi;end

if ty=='lp',
a=sin((fc-f)*pi/sf)/sin((fc+f)*pi/sf);a=a*(abs(a)>10*eps);
na=[1 -a];da=[-a 1];

elseif ty=='hp',
a=-cos((fc+f)*pi/sf)/cos((fc-f)*pi/sf);a=a*(abs(a)>10*eps);
if abs(a)>1,a=1/a;end,
na=[-1 -a];da=[a 1];

elseif ty=='bp',
a=cos(ws)/cos(wd);a=a*(abs(a)>10*eps);k=tan(pi*fc/sf)/tan(wd);
k=k*(abs(k)>10*eps);a1=-2*a*k/(k+1);a2=(k-1)/(k+1);
na=-[1 a1 a2];da=[a2 a1 1];

elseif ty=='bs',
a=cos(ws)/cos(wd);a=a*(abs(a)>10*eps);k=tan(pi*fc/sf)*tan(wd);
k=k*(abs(k)>10*eps);a1=-2*a/(k+1);a2=(1-k)/(k+1);
na=[1 a1 a2];da=[a2 a1 1];

else,error('Unknown filter type'),return,end

else
%wc = 2*pi*fc;
%if wc ~=1
%[n,d]=polymapn(n,d,[wc 0],1);
%end

if ty=='bp'|ty=='bs'
if nargin<7,
a=cos(ws)/cos(wd);else, 
w0=2*pi*fc/sf;
a=cos(w0);
end,end


if ty=='lp',
k=1/tan(f*pi/sf);
na=[k -k];da=[1 1];

elseif ty=='hp',
k=tan(f*pi/sf);
na=[k k];da=[1 -1];

elseif ty=='bp',
%a=cos(ws)/cos(wd);
k=1/tan(wd);a=a*(abs(a)>10*eps);k=k*(abs(k)>10*eps);
na=[k -2*a*k  k];da=[1 0 -1];

elseif ty=='bs',
%a=cos(ws)/cos(wd);
k=tan(wd);a=a*(abs(a)>10*eps);k=k*(abs(k)>10*eps);
na=[k 0 -k];da=[1 -2*a 1];

else,error('Unknown filter type'),return,end
end

[nn,dd]=polymap(n,d,na,da);
dd1=dd;while dd1(1)==0,dd1(1)=[];end
nn=nn/dd1(1);dd=dd/dd1(1);