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

pi=3.14159265358979;

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=asin(imag(zd1)/Absd);//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',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=asin(imag(zd1)/Absd);//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);


