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

matlab之经典数字滤波函数介绍
南京理工大学仪器科学与技术专业 谭彩铭
2010-3-12
1 butter函数
设计一个9阶高通Butterworth数字滤波器，截止频率为300Hz
Fs=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 filter
[bz,az]=impinvar(b,a,Fs)

图5

3 bilinear函数
用双线性变换法数字仿真模拟Butterworth滤波器，程序如下。
Fs=4000;
[b,a]=butter(2,2*pi*1e3,'s')%design analog Butterworth lowpass filter
[bz,az]=bilinear(b,a,Fs)

图6

4 fir1函数
由理想滤波器幅频特性反推滤波系数，得出来的系数数量是无穷多的。故可采用加窗的方法舍去部分，留下有限的滤波系数数量，使仍能基本达到需要的滤波效果。
Window=boxcar(8);
b=fir1(7,0.4,Window)
freqz(b,1)

图7

图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')

图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为采样率）。
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 generater
fs=100;
Ts=1/fs;
sampletime=3;
t=0:Ts:sampletime;
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);
%print effect
figure(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 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)')

下面研究下，它的零相位处理是如何做到的？
直观的思维是先执行filter函数，在做平移即可。那么matalb中filtfilt函数是如何实现的呢？
此滤波是对信号x做前向和反向处理。

9 fftfilt函数
此函数基于FFT和重叠相加法的FIR滤波。

9