自编IIR型50hz陷波器,Matlab源代码,有详细注释

Matlab源代码:%%巴特沃斯型50hz陷波器，可改变陷波带宽和阶数function[Num,Den]=ZB_50_filter(f0,B1,N)%输入参数%f0陷波器中心频率%B1单边带宽%N滤波器阶数fs=1000;%采样率T=1/fs;rp=3;%通带衰减rs=N/2*10;%阻带衰减wp1=((f0-B1)/1000)*2*pi;%下通带截止频率wp2=((f0+B1)/1000)*2*pi;%上通带截止频率ws1=((f0-B1/4)/1000)*2*pi;%阻带下限频率ws2=((f0+B1/4)/1000)*2*pi;%阻带上限频率wc1=(2/T)*tan(wp1/2);wc2=(2/T)*tan(wp2/2);wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);w0=sqrt(wc2*wc1);%阻带中心频率B=wc2-wc1;%带宽wp=1;ws=wp*(wr1*B)/(w0^2-wr1^2);%归一化阻带截止频率[N,wc]=buttord(wp,ws,rp,rs,'s');[Z,P,K]=buttap(N);[Bd,Ad]=zp2tf(Z,P,K);%将零极点转化成传输形式[B,A]=lp2bs(Bd,Ad,w0,B);%对低通滤波器进行频率转换，为带阻滤波器[b,a]=bilinear(B,A,fs);%利用双线性变换法，转换成数字滤波器Num=b;%滤波器分子Den=a;%滤波器分母end举例：设计中心频率为50hz，带宽为0.6hz的IIR的陷波器：f0=50;B1=0.3;%单边带宽0.3hzN=2;[Num,Den]=ZB_50_filter(f0,B1,N);%调用函数计算结果：Num=[0.998118588556613-1.898537748170150.998118588556613];Den=[1-1.898537748170150.996237177113225]调用matlab中的fdatool滤波器工具箱验证：与fdatool自带IIR的butterworth滤波器的幅频曲线一致，滤波器系数几乎一致：