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);