filter22.m

function [num,den,om_n,xi_n,om_d,xi_d]=filter22(Mt,ft,xid,Ts);
%function [num,den,om_n,xi_n,om_d,xi_d]=filter22(Mt,ft,xid,Ts);
%digital filter design function. Design bandstop digital filter
%with 2poles and 2zeros
%inputs:
%Mt ... desired attenuation (Mt<0)/amplification (Mt>0) in dB
%ft ... bandstop frequency in Hz
%xid ... damping of denominator poles for original analog filter. 
%        The resulting denominator damping of digital filter will
%        be in general close to this value. Recommended values 
%        are from 0.6 to 0.9
%Ts ... sampling time
%outputs:
%num ... numerator of digital filter
%den ... denominator of digital filter
%om_n ... natural frequency of numerator [rad/s]
%xi_n ... damping of numerator
%om_d ... natural frequency of denominator [rad/s]
%xi_d ... damping of denominator
%
%written by:  H. Prochazka, I.D. Landau
%7th june 2002

om_zt=2*pi*Ts*ft;

om_st=2*tan(om_zt/2)/Ts;

Mlin=10^(Mt/20);

xin=xid*Mlin;

bz0=4/Ts^2 + 4*xin*om_st/Ts + om_st^2;
bz1=2*om_st^2 - 8/Ts^2;
bz2=4/Ts^2 - 4*xin*om_st/Ts + om_st^2;

az0=4/Ts^2 + 4*xid*om_st/Ts + om_st^2;
az1=2*om_st^2 - 8/Ts^2;
az2=4/Ts^2 - 4*xid*om_st/Ts + om_st^2;

num=[bz0 bz1 bz2];
den=[az0 az1 az2];
%polynomials to frequency and damping
%zeros
zd1 = -num(2) + sqrt( num(2)^2-4*num(1)*num(3) );
zd1=zd1/( 2*num(1) );
Absd=abs(zd1);
Phsd=angle(zd1);
om_n=sqrt( (Phsd^2+log(Absd)^2)/Ts^2 );
xi_n=-log(Absd)/(om_n*Ts);
disp('ZEROS:');
disp(sprintf('nat.frequency[rad/s]\t nat.frequency[Hz]\t damping'));
disp(sprintf('%f\t\t %f\t\t %f\n',om_n,om_n/(2*pi),xi_n));
%poles
zd1 = -den(2) + sqrt( den(2)^2-4*den(1)*den(3) );
zd1=zd1/( 2*den(1) );
Absd=abs(zd1);
Phsd=angle(zd1);
om_d=sqrt( (Phsd^2+log(Absd)^2)/Ts^2 );
xi_d=-log(Absd)/(om_d*Ts);
disp('POLES:');
disp(sprintf('nat.frequency[rad/s]\t nat.frequency[Hz]\t damping'));
disp(sprintf('%f\t\t %f\t\t %f',om_d,om_d/(2*pi),xi_d));
%normalization
num=num/den(1);
den=den/den(1);