|
楼主 |
发表于 2010-5-31
|
(三)、滤波器设计:
1、相关原理:
设计数字滤波器的任务就是寻求一个因果稳定的线性时不变系统,并使系统函数H(z)具有指定的频率特性。
数字滤波器从实现的网络结构或者从单位冲激响应分类,可以分成无限长单位冲激响应(IIR)数字滤波器和有限长单位冲激响应(FIR)数字滤波器。
数字滤波器频率响应的三个参数:
(1) 幅度平方响应:
(2) 相位响应
其中,相位响应
(3) 群时延响应
IIR数字滤波器:
IIR数字滤波器的系统函数为 的有理分数,即
IIR数字滤波器的逼近问题就是求解滤波器的系数 和,使得在规定的物理意义上逼近所要求的特性的问题。如果是在s平面上逼近,就得到模拟滤波器,如果是在z平面上逼近,则得到数字滤波器。
FIR数字滤波器:
设FIR的单位脉冲响应h(n)为实数,长度为N,则其z变换和频率响应分别为
按频域采样定理FIR数字滤波器的传输函数H(z)和单位脉冲响应h(n)可由它的N个频域采样值H(k)唯一确定。
MATLAB中提供了几个函数,分别用于实现IIR滤波器和FIR滤波器。
(1)卷积函数conv
卷积函数conv的调用格式为 c=conv(a,b)
该格式可以计算两向量a和b的卷积,可以直接用于对有限长信号采用FIR滤波器的滤波。
(2)函数filter
函数filter的调用格式为 y=filter(b,a,x)
该格式采用数字滤波器对数据进行滤波,既可以用于IIR滤波器,也可以用于FIR滤波器。其中向量b和a分别表示系统函数的分子、分母多项式的系数,若a=1,此时表示FIR滤波器,否则就是IIR滤波器。该函数是利用给出的向量b和a,对x中的数据进行滤波,结果放入向量y。
(3)函数fftfilt
函数fftfilt的调用格式为 y=fftfilt(b,x)
该格式是利用基于FFT的重叠相加法对数据进行滤波,这种频域滤波技术只对FIR滤波器有效。该函数是通过向量b描述的滤波器对x数据进行滤波。
关于用butter函数求系统函数分子与分母系数的几种形式。
[b,a]=butter(N,wc,'high'):设计N阶高通滤波器,wc为它的3dB边缘频率,以 为单位,故 。
[b,a]=butter(N,wc):当wc为具有两个元素的矢量wc=[w1,w2]时,它设计2N阶带通滤波器,3dB通带为,w的单位为 。
[b,a]=butter(N,wc,'stop'):若wc=[w1,w2],则它设计2N阶带阻滤波器,3dB通带为 ,w的单位为。
如果在这个函数输入变元的最后,加一个变元“s”,表示设计的是模拟滤波器。这里不作讨论。
为了设计任意的选项巴特沃斯滤波器,必须知道阶数N和3dB边缘频率矢量wc。这可以直接利用信号处理工具箱中的buttord函数来计算。如果已知滤波器指标, , 和 ,则调用格式为
[N,wc]=buttord(wp,ws,Rp,As)
对于不同类型的滤波器,参数wp和ws有一些限制:对于低通滤波器,wpws;对于带通滤波器,wp和ws分别为具有两个元素的矢量,wp=[wp1,wp2]和ws=[ws1,ws2],并且 ws1
2、设计内容:
(1)滤波器示例:
在这里为了说明如何用MATLAB来实现滤波,特举出一个简单的函数信号滤波实例(对信号x(n)=sin( n/4)+5cos( n/2)进行滤波,信号长度为500点),从中了解滤波的实现过程。程序如下:
Wn=0.2*pi;
N=5;
[b,a]=butter(N,Wn/pi);
n=0:499;
x=sin(pi*n/4)+5*cos(pi*n/2);
X=fft(x,4096);
subplot(221);plot(x);title('滤波前信号的波形');
subplot(222);plot(X);title('滤波前信号的频谱');
y=filter(b,a,x);
Y=fft(y,4096);
subplot(223);plot(y);title('滤波后信号的波形');
subplot(224);plot(Y);title('滤波后信号的频谱');
在这里,是采用了butter命令,设计出一个巴特沃斯低通滤波器,从频谱图中可以很明显的看出来。下面,也就是本课题的主要内容,也都是运用到了butter函数,以便容易的得到系统函数的分子与分母系数,最终以此来实现信号的滤波。
(2)N阶高通滤波器的设计(在这里,以5阶为例,其中wc为其3dB边缘频率,以 为单位),程序设计如下:
x=wavread('ding.wav');
sound(x);
N=5;wc=0.3;
[b,a]=butter(N,wc,'high');
X=fft(x);
subplot(321);plot(x);title('滤波前信号的波形');
subplot(322);plot(X);title('滤波前信号的频谱');
y=filter(b,a,x);
Y=fft(y);
subplot(323);plot(y);title('IIR滤波后信号的波形');
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
z=fftfilt(b,x);
Z=fft(z);
subplot(325);plot(z);title('FIR滤波后信号的波形');
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
得到结果如图:
(3)N阶低通滤波器的设计(在这里,同样以5阶为例,其中wc为其3dB边缘频率,以 为单位),程序设计如下:
x=wavread('ding.wav');
sound(x);
N=5;wc=0.3;
[b,a]=butter(N,wc);
X=fft(x);
subplot(321);plot(x);title('滤波前信号的波形');
subplot(322);plot(X);title('滤波前信号的频谱');
y=filter(b,a,x);
Y=fft(y);
subplot(323);plot(y);title('IIR滤波后信号的波形');
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
z=fftfilt(b,x);
Z=fft(z);
subplot(325);plot(z);title('FIR滤波后信号的波形');
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
得到结果如图:
(4)2N阶带通滤波器的设计(在这里,以10阶为例,其中wc为其3dB边缘频率,以 为单位,wc=[w1,w2],w1 wc w2),程序设计如下:
x=wavread('ding.wav');
sound(x);
N=5;wc=[0.3,0.6];
[b,a]=butter(N,wc);
X=fft(x);
subplot(321);plot(x);title('滤波前信号的波形');
subplot(322);plot(X);title('滤波前信号的频谱');
y=filter(b,a,x);
Y=fft(y);
subplot(323);plot(y);title('IIR滤波后信号的波形');
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
z=fftfilt(b,x);
Z=fft(z);
subplot(325);plot(z);title('FIR滤波后信号的波形');
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
得到结果如图:
(5)2N阶带阻滤波器的设计(在这里,以10阶为例,其中wc为其3dB边缘频率,以 为单位,wc=[w1,w2],w1 wc w2),程序设计如下:
x=wavread('ding.wav');
sound(x);
N=5;wc=[0.2,0.7];
[b,a]=butter(N,wc,'stop');
X=fft(x);
subplot(321);plot(x);title('滤波前信号的波形');
subplot(322);plot(X);title('滤波前信号的频谱');
y=filter(b,a,x);
Y=fft(y);
subplot(323);plot(y);title('IIR滤波后信号的波形');
subplot(324);plot(Y);title('IIR滤波后信号的频谱');
z=fftfilt(b,x);
Z=fft(z);
subplot(325);plot(z);title('FIR滤波后信号的波形');
subplot(326);plot(Z);title('FIR滤波后信号的频谱');
得到结果如图: |
|