function [na,da] = blt2ord(ty,sf,df,f0)
% BLT2ORD Design 2nd Order peaking or notch IIR digital filters.
%
%       [N,D] = BLT2ORD(TY,SF,DF,F0) Design 2nd ORDER IIR peaking or notch  
%	        filters using the bilinear transformation.
%	TY = 'bp'(bandpass or peaking), or 'bs'(bandstop or notch),  
%	SF = sampling frequency in HERTZ
%	If SF has two elements, the second equals the attenuation A (dB). 
%	DF = A-dB Bandwidth (Hz)   F0=center frequency (Hz) 
%	If DF=[F1 F2] is a 2-element array of band edges, then F0 is ignored
%	N, D contain the coefficients of H(z) in descending order.
%
%	LP2IIR (with no input arguments) invokes the following example:
%
%	Design a 2nd order peaking filter with sampling frequency S = 500 Hz, 
%       2 dB bandwidth = 5 Hz and center frequency  100 Hz and plot results
%        >>[nb,db] = blt2ord('bp',[600 2],5,100); 
%        >>tfplot('z',nb,db,[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: akambard@mtu.edu
% Copyright (c) 1998



if nargin==0,help blt2ord,disp('Strike a key to see results of last example')
pause,[nb,db] = blt2ord('bp',[500 2],5,100)
tfplot('z',nb,db,[0 0.5],0,1); return,end

if length(sf)==2
e=1/sqrt(10^(0.1*sf(2))-1);sf=sf(1);
else
e=1;
end
if length(df)==2,
ws=sum(df)*pi/sf;wd=abs(diff(df))*pi/sf;
a=cos(ws)/cos(wd);
else, 
wd=pi*df/sf;w0=2*pi*f0/sf;a=cos(w0);
end
C=e*tan(wd);a=a*(abs(a)>10*eps);%k=k*(abs(k)>10*eps);
if ty=='bp',
K=C/(1+C);na=[K 0 -K];da=[1 -2*a/(1+C) (1-C)/(1+C)];
elseif ty=='bs',
K=1/(1+C);na=K*[1 -2*a  1];da=[1 -2*a/(1+C) (1-C)/(1+C)];
else,
error('Unknown filter type'),
return,
end