Matlab设计FIR数字滤波器

FIR
数
字
滤
波
器
专业：
学号：
XX：
一课题目的：
1学会使用Matlab的各项功能。

2学会把自己在课堂上学习的知识运用到实践当中。

3了解利用Matlab设计FIR数字滤波器的基本方法。

4在课程设计的过程中掌握程序编译及软件设计的基本方法。

5提高自己对于新知识的学习能力及进行实际操作的能力。

二课题要求：
在信号处理过程中所处理的信号往往混有噪音，从接受到的信号中消除或减弱噪音是信号处理过程中十分重要的问题。

根据有用信号和噪音的不同特性，提取有用信号的过程称为滤波，实现滤波功能的系统称为滤波器。

而数字滤波器又是滤波器中运用极为广泛的一种滤波器。

数值滤波技术是数字信号处理的一个重要组成部分，滤波器的设计是信号处理的核心问题之一。

FIR数字滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。

要求通过网络及各种资料解决实际问题设计一个符合要求的FIR数字滤波器。

三课题内容：
数字滤波器和模拟滤波器有着相同的滤波概念,根据其频率响应特性可分为低通、高通、带通、带阻等类型。

与模拟滤波器相比,数字滤波器除了具有数字信号处理固有优点外,还有滤波精度高、稳定性好、灵活性强等优点。

在数字信号处理中,由于信号中经常混有各种复杂成分，所以很多信号分析都是基于滤波器而进行的，FIR数字滤波器在数字信号处理中发挥着重要作用，采用Matlab软件对FIR数字滤波器进行仿真设计，简化了设计中繁琐的计算。

设计中采用窗函数法，频率采样法和优化设计方法，通过调用Matlab函数设计FIR数字滤波器。

绘制出滤波器的特性图。

利用所设计的滤波器对多个频带叠加的正弦信号进行处理，对比滤波前后的信号时域和频域图，验证滤波器的效果。

最后录制一段语音信号，并对录制的信号进行采样和加噪，绘制出采样后语音信号的时域波形和频谱图，然后用所设计的滤波器对加噪后的信号进行滤波，绘制出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化。

1）窗函数法的Matlab实现
设计FIR数字滤波器的最简单的方法就是窗函数法，通常也称之为傅里叶级数法，FIR数字滤波器的设计首先给出要求的理想滤波器的频率响应H d(e jw),设计一个FIR数字滤波器频率响应H d(e jw)，去逼近理想的滤波响应H d(e jw)。

然而窗函数法设计FIR数字滤波器是在时域进行的，因而必须由理想的频率响应H d(e jw)推导出对应的单位取样响(n)，再设计一个FIR数字滤波器的单位取样响应h(n)去逼近应h
d
(n)。

h
d
窗函数主要用来减少序列因截断而产生的Gibbs效应，但当这个窗函数为矩形时，得到的FIR数字滤波器幅频响应会有明显的Gibbs
效应，并且任意的增加窗函数的长度，Gibbs效应也不能得到改善。

为了克服这种现象，窗函数应该使设计的滤波器具有以下几点：（1）频率特性的主瓣宽度应该尽量窄，且尽可能地将能量集中在主瓣内；
（2）窗函数频率特性的旁瓣在旁瓣ω趋于π的过程中，其能量迅速减小为零。

程序中fir1函数的用法：b=fir1(n,Wn,’ftype’,window)
①n为滤波器的阶数
②Wn为滤波器的截止频率，它是一个0到1的数。

如果Wn是
一个含有两个数的向量，则函数返回一个带通滤波器
③ftype为滤波器的类型，ftype=’high’时，设计的是高通滤波
器；ftype=’stop’时，设计的是带阻滤波器；没有此参数时，
设计的是低通滤波器
④window为指定的窗函数，矩形窗为boxcar(n)，汉宁窗为
hanning(n)，海明窗为hamming(n)，布莱克曼窗为
blackman(n)，凯撒窗为kaiser(n,beta)，没有此参数时，默认为
hamming窗
函数程序如下：
f1=100;
f2=200; %待滤波正弦信号频率
fs=2000; %采样频率
m=(0.3*f1)/(fs/2); %定义过度带宽
M=round(8/m); %定义窗函数的长度
N=M-1; %定义滤波器的阶数
b=fir1(N,0.5*f2/(fs/2)); %使用fir1函数设计滤波器%输入的参数分别是滤波器的阶数和截止频率
figure(1)
[h,f]=freqz(b,1,512); %滤波器的幅频特性图
plot(f*fs/(2*pi),20*log10(abs(h))) %参数分别是频率与幅值
xlabel('频率/赫兹');
ylabel('增益/分贝');
title('滤波器的增益响应');
figure(2)
subplot(211)
t=0:1/fs:0.5; %定义时间X围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %滤波前信号
plot(t,s); %滤波前的信号图像
xlabel('时间/秒');
ylabel('幅度');
title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512); %将信号变换到频域
AFs=abs(Fs); %信号频域图的幅值
f=(0:255)*fs/512; %频率采样
plot(f,AFs(1:256)); %滤波前的信号频域图
xlabel('频率/赫兹');
ylabel('幅度');
title('信号滤波前频域图');
figure(3)
sf=filter(b,1,s);%使用filter函数对信号进行滤波subplot(211)
plot(t,sf)%滤波后的信号图像
xlabel('时间/秒');
ylabel('幅度');
title('信号滤波后时域图');
axis([0.20.5-22]);%限定图像坐标X围subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图xlabel('频率/赫兹');
ylabel('幅度');
title('信号滤波后频域图');
2)频率抽样法的Matlab实现
频率采样法是从频域出发，根据频域采样定理的频率响应H d(e jw)加以
(k)：
等间距抽样，得到h
d
H d(k)=H d(e jw)|ω=(2π)k/N(k=0,1,2...N-1)
(k)可求得FIR滤波器的系统函数H(Z)及频率响应H d(e jw)。

再利用H
d
而在各采样点间的频率响应则是其的加权内插函数延伸叠加的结
果。

但对于一个无限长的序列，用频率采样法必然有一定的逼近误
差，误差的大小取决于理想频响曲线的形状，理想频响特性变换越平缓，则内插函数值越接近理想值，误差越小。

为了提高逼近的质量，可以通过在频率相应的过度内插入比较连续的采样点，扩展过渡带使其比较连续，从而使得通带和阻带之间变换比较缓慢，以达到减少逼近误差的目的。

增大阻带衰减的三种方法：
a)加宽过渡带带宽，以牺牲过渡带换取阻带衰减的增加。

b)过渡带的优化设计。

c)增大N。

直接从频域进行设计，物理概念清楚，直观方便；适合于窄带滤波器的设计，这时频率响应只有少数几个非零值，但是截止频率难以控制。

函数程序如下：
wp=0.2*pi;
wr=0.4*pi;
tr_width=wr-wp;
N=ceil(6.6*pi/tr_width)+1
n=0:1:N;
wc=(wr+wp)/2;
hd=ideal_lp(wc,N);
w_ham=(hamming(N))’;
h=hd.*w_ham;
[db,mag,pha,w]=freqz_m(h,[1]);
delta_w=2*pi/1000;
Ap=-(min(db(1:1:wp/delta_w+1)))
Ar=-round(max(db(wr/delta_w+1:1:501)))
n(:,35)=[];
subplot(2,2,1);stem(n,hd);title(‘理想单位脉冲响应’)
subplot(2,2,2);stem(n,w_ham);title(‘海明窗’)
subplot(2,2,3);stem(n,h);title(‘实际单位脉冲响应’)
subplot(2,2,4);stem(w/pi,db);title(‘幅度响应’)
axis([0,1,-100,10]);
3)最优化设计的Matlab实现
最优化设计方法是指采用最优化准则来设计的方法。

在FIR DF 的最优化设计中，最优化准则有均方误差最小化准则和等波纹切比
雪夫逼近准则两种。

实际设计中，只有采用窗函数才能满足前一种
最优化准则，但由于Gibbs效应的存在，使其根本不能满足设计的要求。

为了满足设计的要求，可以采用其他的窗函数来消除Gibbs效
应，但此时的设计已经能满足该最优化准则了。

因此，要完成FIR
DF的最优化设计，只能采用后一种优化准则来实现。

尽管窗函数法与频率采样法在FIR数字滤波器的设计中有着广泛的应用，但两者不是最优化的设计。

通常线性相位滤波在不同的频
带内逼近的最大容许误差要求不同。

等波纹切比雪夫逼近准则就是
通过通带和阻带使用不同的加权函数，实现在不同频带的加权误差
最大值相同，从而实现其最大误差满足性能指标的条件下达到最大
值，即使得H d(e jw)和H(e jw)之间的绝对误差最小。

尽管按照FIR数字滤波器单位取样响应h(n)的对称性和N的奇，偶性，FIR数字滤波器可以分为4种类型，但是滤波器的频率响应可
以写成统一的形式：
H(e jw)= e-j（N-1）w/2 e j(π/2)k H(ω)
其中，k∈{0,1},H(ω)为幅度函数，且是一个纯实数，表达式也可以写
成统一的形式：
H d(e jw)=Q(ω)P(ω)
其中，Q(ω)为ω的固定函数，P(ω)为M个余弦函数的线性组合。

在优化设计的Matlab实现中，程序中经常使用remez函数，这种函数的使用方法为：
b=remez(n,f,a,w,’ftype’)
1)n为待设计滤波器的阶数；f是一个向量，它是一个0到1的
正数；
2)a是一个向量，指定频率段的幅度值；w对应于各个频段的加
权值；
3)函数的返回值b是设计出的滤波器的系数组成的一个长度为
n+1的向量。

利用Remez函数设计等波纹低通滤波器
设计要求：
1) 通带截频0.5，阻带截频0.6,采样频率2000Hz
2) 带衰减大于等于40dB，通带波纹0.1710和阻带波纹0.01
函数程序如下：
fs=2000; %设定采样频率
rp=3; %通带波纹
rs=40;%阻带波纹
f=[500600];%截止频率
a=[10]; %期望幅度
dev=[(10^(rp/20)-1)/(10^(rp/20)+1)10^(-rs/20)];
[n,fo,ao,w]=remezord(f,a,dev,fs);
b=remez(n,fo,ao,w);
figure(1)
freqz(b,1,1024,fs); %滤波器的特性图
f1=400;
f2=700 ; %待滤波正弦信号频率
t=0:1/fs:0.1; %定义时间X围和步长
s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %滤波前信号
figure(2)
subplot(211)
plot(t,s); %滤波前的信号图像
xlabel('时间/秒');
ylabel('幅度');
title('信号滤波前时域图');
subplot(212)
Fs=fft(s,512); %将信号变换到频域
AFs=abs(Fs); %信号频域图的幅值
f=(0:255)*fs/512; %频率采样
plot(f,AFs(1:256)); %滤波前的信号频域图
xlabel('频率/赫兹');
ylabel('幅度');
title('信号滤波前频域图');
figure(3)
sf=filter(b,1,s); %使用filter函数对信号进行滤波subplot(211)
plot(t,sf) %滤波后的信号图像
xlabel('时间/秒');
ylabel('幅度');
title('信号滤波后时域图');
subplot(212)
Fsf=fft(sf,512);%滤波后的信号频域图
AFsf=abs(Fsf);%信号频域图的幅值
f=(0:255)*fs/512;%频率采样
plot(f,AFsf(1:256))%滤波后的信号频域图
xlabel('频率/赫兹');
ylabel('幅度');
title('信号滤波后频域图');
四实验心得：
1通过这个让我更灵活的运用Matlab来设计，也让我学会了FIR的各种设计方法，我相信这会是我人生很重要的一课，使我受益良多，提高自己对于新知识的学习能力及进行实际操作的能力。

让我学会了在不懂的时候也不要放弃，要去找方法来解决问题，才能够收获知识，让我们更进一步的成长，学会在问题中找方法，找出问题的症结，才能解决问题，不然会一直的困在问题中不能前进，要相信自己能行。

这次试验中我受益良多，学到了很多东西。