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

f=5;
x=sin(2*pi*f*t);
%filter generater
fc=4;
%cutoff frequency in now sample rate
b=fir1(40,fc/(0.5*fs));
figure(1);
freqz(b,1,128,fs)
%filter course
y=filter(b,1,x);
Fs=4000; [b,a]=butter(2,2*pi*1e3,'s')%design analog Butterworth lowpass filter [bz,az]=bilinear(b,a,Fs)
图6 4 fir1 函数 由理想滤波器幅频特性反推滤波系数，得出来的系数数量是无穷多的。故可采用加窗的方法 舍去部分，留下有限的滤波系数数量，使仍能基本达到需要的滤波效果。
matlab 之经典数字滤波函数介绍 南京理工大学仪器科学与技术专业 谭彩铭
2010-3-12 1 butter 函数 设计一个 9 阶高通 Butterworth 数字滤波器，截止频率为 300Hz
Fs=1000; [b,a]=butter(9,300/(Fs/2),'high') freqz(b,a,128,Fs)
fs=100; fc=[0.35 0.45]; b=fir1(40,fc); figure(1); freqz(b,1,128,fs)
图 11
5
图 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)。如下图程序所示。
Window=boxcar(8); b=fir1(7,0.4,Window) freqz(b,1)
图7
3
图8 5 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')
图3
图4 可见 h1 和 h 相等 图 2 中出现小于-360 度是否表达其他不同的意义？-361 度和-1 度有什么区别吗？对于正弦 波应该是一样的，故理论上说应该没有区别。 butter 函数的原理是什么？ 顾名思义，butter 函数的原理是基于 Butterworth 滤波器。 这里始终要带着这个问题去研究，滤波系数本身有什么特性竟然可使低频的滤掉，高频的通 过，其实这里想要寻找的是敏捷控程，理论上这个问题的答案已经很成熟。 2 impinvar 函数 用冲击响应不变法数字仿真模拟 Butterworth 滤波器，程序如下。
ቤተ መጻሕፍቲ ባይዱ图1
图2 下面看一看 freqz 函数
例如对离散系统传递函数
H
(z)

b(1) a(1)

b(2)z 1 a(2)z 1
b(3)z 2 a(3)z 2
freqz 函数的的主要计算环节是计算 H (e jwTs )
编写下图所示程序验证之
1
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)')
下面研究下，它的零相位处理是如何做到的？ 直观的思维是先执行 filter 函数，在做平移即可。那么 matalb 中 filtfilt 函数是如何实现的呢？ 此滤波是对信号 x 做前向和反向处理。 9 fftfilt 函数 此函数基于 FFT 和重叠相加法的 FIR 滤波。
9
- a(2)*y(n-1) - ... - a(na+1)*y(n-na) 8 filtfilt 函数
8
Zero-phase digital filtering filtfilt 函数的计算方式同 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(x)','filtfilt(x)')
2
Fs=4000; [b,a]=butter(2,2*pi*1e3,'s') %design analog Butterworth lowpass filter [bz,az]=impinvar(b,a,Fs)
图5 3 bilinear 函数 用双线性变换法数字仿真模拟 Butterworth 滤波器，程序如下。
图9
4
图 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 为采样率）。
fs=100; fc=[20 25]/(0.5*fs); b=fir1(40,fc); figure(1); freqz(b,1,128,fs)
图 13
6
下列程序是较为综合的一个程序
图 14
clear;
%signal generater
fs=100;
Ts=1/fs;
sampletime=3;
t=0:Ts:sampletime;
%print effect
figure(2);
plot(t,x);
hold on;
plot(t,y,'r');
hold off;
7
图 15
图 16
图 17
由图 17 中，可得到信号频率为 5Hz 时的衰减系数约为 0.25， 20log10 (0.25) 12 ，和图
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)