matlab之经典数字滤波函数介绍（副本）.docVIP

matlab之经典数字滤波函数介绍 南京理工大学仪器科学与技术专业 谭彩铭 2010-3-12 1 butter函数 设计一个9阶高通Butterworth数字滤波器，截止频率为300Hz 图1 图2 下面看一看freqz函数 例如对离散系统传递函数 freqz函数的的主要计算环节是计算 编写下图所示程序验证之 图3 图4 可见h1和h相等 图2中出现小于-360度是否表达其他不同的意义？-361度和-1度有什么区别吗？对于正弦波应该是一样的，故理论上说应该没有区别。 butter函数的原理是什么？ 顾名思义，butter函数的原理是基于Butterworth滤波器。 这里始终要带着这个问题去研究，滤波系数本身有什么特性竟然可使低频的滤掉，高频的通过，其实这里想要寻找的是敏捷控程，理论上这个问题的答案已经很成熟。 2 impinvar函数 用冲击响应不变法数字仿真模拟Butterworth滤波器，程序如下。 图5 3 bilinear函数 用双线性变换法数字仿真模拟Butterworth滤波器，程序如下。 图6 4 fir1函数 由理想滤波器幅频特性反推滤波系数，得出来的系数数量是无穷多的。故可采用加窗的方法舍去部分，留下有限的滤波系数数量，使仍能基本达到需要的滤波效果。 图7 图8 5 fir2函数 fir2函数的基本原理同fir1函数，它的功能更进一层，可以设计任意形状的频率响应图形。 图9 图10、 图9中，w如果是角频率值，将w转换成频率值时，应该是将w除以2*pi，但是程序中为什么除以的是pi呢？ 准确地说，图9中的w并非是角频率值，而是频率值，freqz函数调用时若没有加入采样率参数，其返回的频率值的范围是0~pi。 6 fir1函数补充 对调用方式b=fir(n,wn)，wn的取值范围是（0，1），其中1对应于0.5fs（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)。如下图程序所示。 图13 图14 下列程序是较为综合的一个程序 图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 filtering filtfilt函数的计算方式同filter函数一样，不同之处在于做了零相位处理。下图所示程序就可以清楚地看到这点区别。 下面研究下，它的零相位处理是如何做到的？ 直观的思维是先执行filter函数，在做平移即可。那么matalb中filtfilt函数是如何实现的呢？ 此滤波是对信号x做前向和反向处理。 9 fftfilt函数 此函数基于FFT和重叠相加法的FIR滤波。 8 Fs=1000; [b,a]=butter(9,300/(Fs/2),high) freqz(b,a,128,Fs) 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)) Fs=4000; [b,a]=butter(2,2*pi*1e3,s) %design analog Butterworth lowpass filter [bz,az]=impinvar(b,a,Fs) Fs=4000; [b,a]=butter(2,2*pi*1e3,s)%design analog Butterworth lowpass filter [bz,az]=