matlab之经典数字滤波函数介绍.doc

matlab之经典数字滤波函数介绍南京理工大学仪器科学与技术专业 谭彩铭2010-3-121 butter函数设计一个9阶高通Butterworth数字滤波器，截止频率为300HzFs=1000;b,a=butter(9,300/(Fs/2),high)freqz(b,a,128,Fs)图1图2下面看一看freqz函数例如对离散系统传递函数freqz函数的的主要计算环节是计算编写下图所示程序验证之Fs=1000;b,a=butter(2,300/(Fs/2),high)h,f=freqz(b,a,128,Fs);h1=(b(1)+b(2)*exp(-j*2*pi*f/Fs)+b(3)*exp(-2*j*2*pi*f/Fs)./(a(1)+a(2)*exp(-j*2*pi*f/Fs)+a(3)*exp(-2*j*2*pi*f/Fs);plot(abs(h1-h)title(abs(h1-h)图3图4可见h1和h相等图2中出现小于-360度是否表达其他不同的意义？-361度和-1度有什么区别吗？对于正弦波应该是一样的，故理论上说应该没有区别。butter函数的原理是什么？顾名思义，butter函数的原理是基于Butterworth滤波器。这里始终要带着这个问题去研究，滤波系数本身有什么特性竟然可使低频的滤掉，高频的通过，其实这里想要寻找的是敏捷控程，理论上这个问题的答案已经很成熟。2 impinvar函数用冲击响应不变法数字仿真模拟Butterworth滤波器，程序如下。Fs=4000;b,a=butter(2,2*pi*1e3,s) %design analog Butterworth lowpass filterbz,az=impinvar(b,a,Fs)图53 bilinear函数用双线性变换法数字仿真模拟Butterworth滤波器，程序如下。Fs=4000;b,a=butter(2,2*pi*1e3,s)%design analog Butterworth lowpass filterbz,az=bilinear(b,a,Fs)图64 fir1函数由理想滤波器幅频特性反推滤波系数，得出来的系数数量是无穷多的。故可采用加窗的方法舍去部分，留下有限的滤波系数数量，使仍能基本达到需要的滤波效果。Window=boxcar(8);b=fir1(7,0.4,Window)freqz(b,1)图7图85 fir2函数fir2函数的基本原理同fir1函数，它的功能更进一层，可以设计任意形状的频率响应图形。f = 0:0.1:1; m = 0 0 1 1 0 0 1 1 1 0 0;b = fir2(30,f,m);h,w = freqz(b,1,128);plot(f,m,w/pi,abs(h)legend(Ideal,fir2 Designed)title(Comparison of Frequency Response Magnitudes)图9图10、图9中，w如果是角频率值，将w转换成频率值时，应该是将w除以2*pi，但是程序中为什么除以的是pi呢？准确地说，图9中的w并非是角频率值，而是频率值，freqz函数调用时若没有加入采样率参数，其返回的频率值的范围是0pi。6 fir1函数补充对调用方式b=fir(n,wn)，wn的取值范围是（0，1），其中1对应于0.5fs（fs为采样率）。fs=100;fc=0.35 0.45;b=fir1(40,fc);figure(1);freqz(b,1,128,fs)图11图12对图11中程序，1对应0.5fs，那么0.45对应0.45*0.5*fs=22.5Hz，0.35对应0.35*0.5*fs=17.5Hz，和图12吻合。但是实际应用时，我们知道采样率，知道需要的通带，我们需要的是直接的答案。假设采样率为100Hz，现在要设计的带通滤波器的带宽是20Hz，25Hz，对应图11中的fc我们该怎么取值，显然应该是fc=20 25/(0.5*fs)。如下图程序所示。fs=100;fc=20 25/(0.5*fs);b=fir1(40,fc);figure(1);freqz(b,1,128,fs)图13图14下列程序是较为综合的一个程序clear;%signal generaterfs=100;Ts=1/fs;sampletime=3;t=0:Ts:sampletime;f=5;x=sin(2*pi*f*t);%filter generaterfc=4; %cutoff frequency in now sample rateb=fir1(40,fc/(0.5*fs);figure(1);freqz(b,1,128,fs)%filter coursey=filter(b,1,x);%print effectfigure(2);plot(t,x);hold on;plot(t,y,r);hold off;图15图16图17由图17中，可得到信号频率为5Hz时的衰减系数约为0.25，和图16基本吻合。7 filter函数Filter data with an infinite impulse response (IIR) or finite impulse response (FIR) filter对y = filter(b,a,x)，y(n) = b(1)*x(n) + b(2)*x(n-1) + . + b(nb+1)*x(n-nb) - a(2)*y(n-1) - . - a(na+1)*y(n-na)8 filtfilt函数Zero-phase digital filteringfiltfilt函数的计算方式同filter函数一样，不同之处在于做了零相位处理。下图所示程序就可以清楚地看到这点区别。clear;t=0:0.001:0.1;x1=sin(2*pi*40*t);x2=0.5*rand(size(t);x=x1+x2;A=1 -1.143 0.4128;B=0.06745 0.1348 0.06745;y=filter(B,A,x);z=filtfilt(B,A,x);plot(t,x,t,y,t,z);legend(x,filter(