数字滤波器的设计及实现
摘要：通过MATLAB软件的信号处理工具箱中的滤波器设计各种IIR和FIR数字滤波器，并绘制滤波器的幅频特性、相频特性。通过观察滤波器的输入、输出信号的时域波形及其频谱，建立了数字滤波的概念。 关键词：MATLAB；IIR数字滤波器；FIR数字滤波器
（一）设计目的
1、熟悉IIR数字滤波器和FIR数字滤波器的设计原理和方法；
2、学会调用MATLAB信号处理工具箱中的滤波器设计函数设计各种IIR和FIR数字滤波器，学会根据滤波要求确定滤波器指标参数；
3、掌握用IIR和FIR数字滤波器的MATLAB实现方法，并能绘制滤波器的幅频特性、相频特性；
4、通过观察滤波器的输入、输出信号的时域波形及其频谱，建立数字滤波的概念。
（二）设计要求
用MATLAB软件设计IIR数字滤波器和FIR数字滤波器，并绘制滤波器的幅频特性、相频特性。
（三）设计原理
数字滤波器是将输入数字序列通过一定的运算后转变为输出数字序列的数字信号处理器。数字滤波器的输入、输出均为数字信号, 通过一定运算系改变输入信号所含频率成分的相对例或消除某些频率成分。与模拟滤波器相比, 数字滤波器的主要优点是：（1）精度和稳定性高; (2)系统函数容易改变, 因而灵活性高; （3）不存在阻抗匹配问题; ( 4)便于大规模集成; ( 5)可以实现多维滤波。它不仅能实现模拟处理的大部分功能, 而且还能完成模拟处理由于成本、可靠性等原因而无法具体实现的功能。
所谓抑制载波单频调制信号，就是两个正弦信号相乘，它有2个频率成分：和频fc+f0，差频fc-f0，这两个频率成分关于载波频率fc对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率fc对称的两根谱线。显然，当调制频率f0和（或）载波频率fc不同时，可以得到包含不同频率成分的单频调幅信号。
（四）设计内容
1、调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，观察st
的时域波形和幅频特性曲线；
根据题目的要求编写代码如下：
function st=mstg %产生信号序列st，并显示st的时域波形和频谱
%st=mstg返回三路调幅信号相加形成的混合信号，长度N=800
1
信号长度N为800
Fs=10000;T=1/Fs;Tp=N*T; %采样频率Fs=10kHz，Tp为采样时间 t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10; %第1路调幅信号载波频率fc1=1000Hz
fm1=fc1/10; %第1路调幅信号的调制信号频率fm1=100Hz fc2=Fs/20; %第2路调幅信号载波频率fc2=500Hz
fm2=fc2/10; %第2路调幅信号的调制信号频率fm2=50Hz fc3=Fs/40; %第3路调幅信号载波频率fc3=250Hz
fm3=fc3/10; %第3路调幅信号的调制信号频率fm3=25Hz xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t); %产生第1路调幅信号
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t); %产生第2路调幅信号
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t); %产生第3路调幅信号
st=xt1+xt2+xt3; %三路信号相加，得到复合信号
fxt=fft(st,N); %计算信号st的频谱
%以下为绘图命令
subplot(2,1,1);
plot(t,st);grid;xlabel('t/s');ylabel('s(t)');
axis([0,Tp,min(st),max(st)]);title('(a)s(t)的波形')
subplot(2,1,2);
stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b)s(t)的频谱')
axis([0,Fs/8,0,1.2]);
xlabel('f/Hz');ylabel('幅度
');
运行的波形如下：
图一三路调幅信号st的时域波形和幅频特性曲线
2
、 通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率；假定要求滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB，编程调用MATLAB滤波器设计函数分别设计这三个数字滤波器，并绘图显示其幅频特性曲线。
根据题目编写代码如下：
% 低通 巴特沃斯
clear;clc
Fs=10000; %采样频率
fp=280; %通带截止频率
fs=525; %阻带截止频率
rp=0.1; %通带最大衰减
rs=60; %阻带最小衰减
wp=2*fp/Fs;ws=2*fs/Fs; %计算数字滤波器的设计指标
[N,wc]=buttord(wp,ws,rp,rs); %计算数字滤波器的阶数和通带截止频率
[b,a]=butter(N,wc); %计算数字滤波器系统函数
w=0:0.01*pi:pi;
[h,w]=freqz(b,a,w); %计算数字滤波器的幅频响应
h=20*log10(abs(h)); %求频率的幅度值
%绘图程序
subplot(3,1,1);plot(w/pi,h);grid;axis([0,1,-700,40]);
xlabel('\omega/\pi');ylabel('幅度/dB');title('巴特沃斯低通滤波器的幅频特性曲线');
%带通 切比雪夫II
clear;clc
Fs=10000; fp1=400;fp2=600;fs1=300;fs2=750; rp=0.1;rs=60;
wp=[2*fp1/Fs,2*fp2/Fs];ws=[2*fs1/Fs,2*fs2/Fs]; %计算数字滤波器的设计指标
[N,wso]=cheb2ord(wp,ws,rp,rs); %计算数字滤波器的阶数和阻带截止频率
[b,a]=cheby2(N,rs,wso); %计算数字滤波器的系统函数
w=0:0.01*pi:pi;
[h,w]=freqz(b,a,w); %计算数字滤波器的幅频响应
h=20*log10(abs(h));
%绘图程序
subplot(3,1,2);plot(w/pi,h);grid;axis([0,1,-100,50]);
xlabel('\omega/\pi');ylabel('幅度/dB');title('切比雪夫II带通滤波器的幅频特性曲线'); %高通 切比雪夫I
3
Fs=10000;fp=800;fs=600; rp=0.1;rs=60;wp=2*fp/Fs;ws=2*fs/Fs;
[N,wpo]=cheb1ord(wp,ws,rp,rs); %计算数字滤波器的阶数和通带截止频率
[b,a]=cheby1(N,rp,wpo,'high'); %计算数字滤波器系统函数
w=0:0.01*pi:pi;
[h,w]=freqz(b,a,w); %计算数字滤波器的幅频响应
%绘图程序
subplot(3,1,3);plot(w/pi,h);grid;axis([0,1,-250,50]);
xlabel('\omega/\pi');ylabel('幅度/dB');title('切比雪夫I高通滤波器的幅频特性曲线
');
运行波形如下：
图二
3、 用所设计的三个滤波器分别对复合信号st进行滤波，分离出st中的三路不同载波频率的
调幅信号，并绘图显示滤波后信号的时域波形和频谱，观察分离效果。
根据题目编写代码如下：
%产生调幅信号（同内容1，这里省略）
%低通滤波器
fp=300;fs=320;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）
[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和A
4
滤波器软件实现
y1=fft(y1t);
% 绘图部分
figure(2)
subplot(311)
[H,W]=freqz(B,A,800);
plot(W*Fs/2/pi,abs(H));
xlabel('Hz');
ylabel('H(w)');
title('低通滤波器');
axis([0,2000,0,1.2]);
grid;
figure(1)
subplot(4,2,3);
plot(t,y1t);
xlabel('t');
ylabel('y(t)');
title('分离出的250Hz的波形');
figure(1)
subplot(4,2,4);
stem(f,abs(y1)/max(abs(y1)),'.');
xlabel('Hz');
ylabel('|H|');
title('250Hz的频谱');
axis([0,1200,0,1]);
grid;
%带通滤波器
fpl=400;fpu=580;fsl=300;fsu=700;
wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;
[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和A
y2t=filter(B,A,st); %滤波器软件实现
y2=fft(y2t);
%绘图部分
figure(2)
subplot(312)
[H,W]=freqz(B,A,800);
plot(W*Fs/2/pi,abs(H));
xlabel('Hz');
ylabel('H(w)');
title('带通滤波器');
axis([0,2000,0,1.2]);
5
figure(1)
subplot(4,2,5);
plot(t,y2t);
xlabel('t');
ylabel('y(t)');
title('分离出的500Hz的波形');
figure(1)
subplot(4,2,6);
stem(f,abs(y2)/max(abs(y2)),'.'); xlabel('Hz');
ylabel('|H|');
title('频谱');
axis([0,1200,0,1]);
grid;
%高通滤波器
fp=800;fs=780;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF
[N,wp]=ellipord(wp,ws,rp,rs); %止频率wp
[B,A]=ellip(N,rp,rs,wp,'high'); %量B和A
y3t=filter(B,A,st); %y3=fft(y3t);
%绘图部分
figure(2)
subplot(313)
[H,W]=freqz(B,A,800);
plot(W*Fs/2/pi,abs(H));
xlabel('Hz');
ylabel('H(w)');
title('高通滤波器');
axis([0,2000,0,1.2]);
grid;
figure(1)
subplot(4,2,7);
plot(t,y3t);
xlabel('t');
ylabel('y(t)');
title('分离出的1000Hz的波形');
figure(1)
subplot(4,2,8);
stem(f,abs(y3)/max(abs(y3)),'.'); xlabel('Hz');
指标（低通滤波器的通、阻带边界频） 调用ellipord计算椭圆DF阶数N和通带截调用ellip计算椭圆带通DF系统函数系数向滤波器软件实现 6
title('频谱');
axis([0,1200,0,1]); grid;
运行波形如下：
图三
（五）调试分析
由图一可见，三路信号时域混叠无法在时域进行分离，但频域是分离的。容易看出，这三路调幅信号的载波频率分别为250Hz、500Hz和1000Hz，因此可以通过设计合适的滤波器的方法在频域分离。
如图三中经过滤波器的方法在频移分离出了250Hz、500Hz和1000Hz的波形。如图三中在250Hz、500Hz和1000Hz的地方就有被分离出来的波形。
（六）心得体会
通过这次的课程设计使我进一步的熟悉了数字滤波器的原理和设计方法及实现方法，以及一些MATLAB的代码的意思和使用。也使我学会自己搜索资料的能力，为以后的学习打下基础。 参考文献：
[ 1] 任志刚.“数字信号处理”多媒体教学方法初探[J].电气电子教学学报,2006,28(6):102-104。
[ 2] 刘大年,史旺旺,孙贵根,等.“数字信号处理”课程的形象化教学方法探索[ J].电气电子报,2006, 28(4):104-108。
[ 3] 曾孟雄,吴海华.基于LabV IEW 平台的虚拟滤波器设计[J].三峡大学学报:自然科版,2003,25(2):53- 156。
[ 4] 杨乐平,李海涛,赵勇,等.LabV IEW 高级程序设计[M].北京:清华大学出版社,2003。
[ 5] 唐建锋.MATLAB在数字滤波器设计教学中的应用研究[J].衡阳师范学院学报,2005,26(3):164- 166。
7
第二篇：课程设计之数字信号处理
课程设计
————张前杰
Wp=0.2*pi;Ws=0.05*pi;Ap=3;As=48; Wp=tan(Wp/2);Ws=tan(Ws/2);
[N,wn]=buttord(Wp,Ws,Ap,As,'s'); fprintf('滤波器阶数N=%.0f\n',N);
[b,a]=butter(N,1,'s');
[numa,dena]=lp2hp(b,a,Wp);
[numd,dend]=bilinear(numa,dena,0.5); disp('分子系数b'); fprintf('%.4e ',numd); fprintf('\n');
disp('分母系数a'); fprintf('%.4e',dend); fprintf('\n');
w=linspace(0,pi,1024); h=freqz(numd,dend,w); plot(w/pi,20*log10(abs(h))); axis([0,1,-50,0]);grid; xlabel('归一化频率'); ylabel('幅度/dB');
