IIR数字滤波器的设计及软件实现

一．实验目的

（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法; （2）学会用MATLAB信号处理工具箱中的滤波器设计函数（或滤波器设计分析工具FDAtool)设计各种滤波器，学会根据滤波需求确定滤波器指标参数; （3）掌握IIR数字滤波器的MATLAB实现方法；

（4）通过观察滤波器输入、输出信号的时域波形及其频谱，建立数字滤波的概念。

二．实验原理

设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性不变法)，应用最广泛的是双线性变换法。基本的设计过程是：

将给定的数字滤波器指标转换成模拟滤波器的指标; 涉及模拟滤波器；

将模拟滤波器的系统函数转换成数字滤波器的系统函数.

MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。本实验的数字滤波器的MATLAB实验是调用MATLAB信号处理工具箱的函数filter对给定的输入信号x(n）进行滤波,得到滤波后的输出信号y(n). 三．实验内容及步骤

1．信号处产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数还会自动回图显示st的时域波形和幅频特性曲线，由后图可见，三路信号时域混叠无法在时域分离。但频域是分离的，所以可通过滤波的方法在频域分离。

2．将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。且滤波器的通带最大衰减为0。1dB，阻带最小衰减为60bB。

提示：抑制载波单频调幅信号的数学表示式为

1 s(t)cos2ftcos2ftcos2fftcos2fft0cc0c02fc0其中,cos2fct称为载波，f为载波频率，cos2cft称为单频调制信号，f0为调制正弦波信号频率，且满足〉

f0。由上式可见，所谓抑制载波单频调

幅信号，就是2个正弦信号相乘，它有2个频率成分：和频

fc+

f0和差频

fc—

f0，这2个频率成分关于载波频率

fc对称。所以，1路抑制载波单频调幅信号

的频谱图是关于载波频率

fc对称的2根谱线，其中没有载频成分,故取名为抑制

载波单频调幅信号。图中三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz.

0

3. 编程调用MATLAB滤波器涉及函数ellipord和ellip分别设计这三个椭圆滤波器，并绘图显示其损耗函数曲线;

4. 调用滤波器实验函数filter，用三个滤波器分别对信号产生函数mstg产生的信号st进行滤波,分离出st中的三路不同载波频率的调幅信号y1(n)、y2(n）、y3(n)。 滤波器参数的选取:

 对于载波频率为250Hz的条幅信号,可以选用低通滤波器分离，其指标为：

通带截止频率阻带截止频率

fp280Hz，通带最大衰减P0.1dB；

ffs450Hz，阻带最小衰减s60dB;

440Hz,

 对于载波频率为500Hz的条幅信号，可以选用带通滤波器分离,其指标为：

通带截止频率阻带截止频率

plfpu560Hz,通带最大衰减p0.1dB；

ffsl275Hz，

fsu900Hz，阻带最小衰减

s60dB；

 对于载波频率为1000Hz的条幅信号，可以选用高通滤波器分离,其指标为：

通带截止频率阻带截止频率

p890Hz，通带最大衰减P0.1dB；

fs550Hz，阻带最小衰减s60dB;

说明：

（1)为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽可能使滤波器过渡带宽尽可能宽；

（2）与信号产生函数mstg相同采用频率Fs=10kHz; （3）为了滤波器阶数最低，选用椭圆滤波器。

四．试验程序框图

调用函数mstg产生st，自动绘图显示st的时域波形和幅频特性曲调用ellipord和ellip分别设计三个椭圆滤波器，并绘图显示其幅频响应特性曲线 调用filter，用三个滤波器分别对信号st进行滤波，分离出三路不同载波频率的调幅信号y1(n),y2(n),和y3(n) 绘图显示y1(n),y2(n)和y3(n)的时域波形 结束 五．思考题及简答

1

1、阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率 答：

第一路调幅信号的载波频率fc1=1000Hz; 第一路调幅信号的调制频率fm1=100Hz； 第二路调幅信号的载波频率fc2=500Hz; 第二路调幅信号的调制频率fm2=50Hz; 第三路调幅信号的载波频率fc3=250Hz； 第三路调幅信号的调制频率fm3=25Hz；

2、信号产生函数mstg中采样点数N=1600，对st进行N点FFT就可以得到6根理想谱线。如果取N=1800，可否得到6根理想谱线？为什么？N=2000呢?请改变采样点数N的值,观察频谱图验证判断是否正确?

答:因为信号st是周期序列,谱分析时要求观察时间为整数倍周期。分析可知，st的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25×400Hz,即在25Hz的正弦波的1个周期中采样400点。所以，当N为400的整数倍时一定为st的整数个周期.因此，采样点数N=1600和N=2000时，对st进行N点FFT可以得到6根理想谱线。如果取N=1800,不是400的整数倍，不能得到6根理想谱线。

（1）N=1600时:

(a)s(t)的波形32s(t)10-100.0020.0040.0060.0080.010.0120.0140.0160.018t/s(b) s(t)的频谱0.021幅度0.5002004006008001000f/HZ12001400160018002000

（2）N=1800时：

2

(a)s(t)的波形32s(t)10-100.0050.01t/s(b) s(t)的频谱10.0150.02幅度0.5002004006008001000f/HZ12001400160018002000（3）N=2000S时：

(a)s(t)的波形32

s(t)10-100.0050.010.0150.020.025t/s(b) s(t)的频谱1幅度0.5002004006008001000f/HZ120014001600180020003、修改信号产生函数mstg，给每路调幅信号加入载波成分，产生调幅（AM）信号，重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别.

AM信号表示式：st

AAcos2ftcos2ft AAdm0cdm

3

取值：Ad10 ,Am5，结果见（附录Ⅰ）

六．实验结果（程序附录Ⅱ）

原信号输出：

(a)s(t)的波形32s(t)10-100.0010.0020.0030.0040.0050.0060.007t/s(b) s(t)的频谱0.0080.0090.011幅度0.5002004006008001000f/HZ12001400160018002000

低通输出：

(a) s(t)损耗函数曲线0-20-40-60-80幅度（dB）00.10.20.30.40.50.6/(b) s(t)的波形0.70.80.9110.50-0.500.010.020.030.040.05t/s(c) s(t)的频谱0.060.07y(n)幅度10.5002004006008001000f/Hz12001400160018002000

带通输出:

4

(a) s(t)损耗函数曲线0-20-40-60-80幅度（dB）00.10.20.30.40.50.6/(b) s(t)的波形0.70.80.9110.50-0.500.010.020.030.040.05t/s(c) s(t)的频谱0.060.07y(n)幅度10.5002004006008001000f/Hz12001400160018002000高通输出：

(a) s(t)损耗函数曲线0-20-40-60-80

幅度（dB）00.10.20.30.40.50.6/(b) s(t)的波形0.70.80.9110.50-0.500.010.020.030.040.05t/s(c) s(t)的频谱0.060.07y(n)幅度10.5002004006008001000f/Hz12001400160018002000

5

附录（Ⅰ)：

原信号输出：

(a)s(t)的波形4020s(t)0-2000.0010.0020.0030.0040.0050.0060.0070.0080.009t/s(b) s(t)的频谱0.011幅度0.5002004006008001000f/HZ12001400160018002000低通输出：

(a) s(t)损耗函数曲线0-20-40-60-80

幅度（dB）00.10.20.30.40.50.6/(b) s(t)的波形0.70.80.91y(n)100-1000.010.020.030.040.05t/s(c) s(t)的频谱0.060.071幅度0.5002004006008001000f/Hz12001400160018002000

6

带通输出：

(a) s(t)损耗函数曲线0-20-40-60-80幅度（dB）00.10.20.30.40.60.5/(b) s(t)的波形0.70.80.91y(n)100-1000.010.020.030.050.04t/s(c) s(t)的频谱0.060.071幅度0.5002004006008001000f/Hz12001400160018002000高通输出:

(a) s(t)损耗函数曲线0-20-40-60-80

幅度（dB）00.10.20.30.40.50.6/(b) s(t)的波形0.70.80.91y(n)100-1000.010.020.030.040.05t/s(c) s(t)的频谱0.060.071幅度0.5002004006008001000f/Hz12001400160018002000

7

附录(Ⅱ）:

%主函数

%IIR数字滤波器设计及软件实现 clear all;

％调用信号产生函数mstg产生又三路抑制载波调幅信号相加构成的复合信号 syms st; syms t；

st=mstg; ％低通滤波器设计与实现 Fs=10000；T=1/Fs;

n=800；Tp=n＊T；k=0：n-1;f=k/Tp； fp=280；fs=450；

wp=2*fp/Fs；ws=2＊fs/Fs;

rp=0.1;rs=60; ％DF指标;(低通滤波器的通阻带边界频率）

［N，wp0]=ellipord（wp，ws,rp，rs);％调用ellipod计算椭圆DF阶数N和通带截止频率wp

［B,A]=ellip(N，rp，rs,wp0）；%调用ellip计算椭圆带通DF系统函数系数向量B和A y1t=filter(B，A，st);％滤波器的软件实现 fyt=fft(y1t，n)； %下面为绘图部分 figure(2）； subplot（3,1,1); myplot（B，A）; yt='y_1(t）’；

subplot（3，1，2)； tplot(y1t,T,yt）; subplot(3,1,3）;

stem(f，abs（fyt）/max（abs（fyt))，'。’）； grid；

title（’（c) s(t)的频谱’); axis([0,Fs/5,0,1.2］）；

xlabel(’f/Hz’)；ylabel（’幅度');

％带通滤波器的实现与设计

fpl=440；fpu=560;fsl=275;fsu=900； wp=[2＊fpl/Fs,2＊fpu/Fs]； ws=[2＊fsl/Fs,2＊fsu/Fs]； rp=0。1；rs=60;

[N，wp0］=ellipord（wp，ws,rp，rs)； ［B，A］=ellip(N,rp，rs，wp0）； y2t=filter（B,A,st）; fyt=fft（y2t,n）； figure(3);

subplot（3,1，1）;

8

myplot(B,A)； yt='y_1(t）’;

subplot（3,1，2）； tplot(y1t,T，yt); subplot(3，1,3)；

stem(f,abs(fyt)/max（abs(fyt)),’.'）; grid；

title('(c) s（t)的频谱’）； axis（［0，Fs/5,0,1。2]);

xlabel（’f/Hz’)；ylabel(’幅度’)；

%高通滤波器的实现与设计 fp=890；fs=600;

wp=2＊fp/Fs;ws=2*fs/Fs; rp=0。1；rs=60;

［N，wp0］=ellipord（wp,ws,rp，rs）； ［B,A］=ellip(N，rp，rs，wp0，’high’）; y3t=filter(B,A，st）; fyt=fft(y3t，n）； figure(4);

subplot(3,1,1）； myplot(B，A）; yt=’y_1(t)';

subplot(3，1，2)； tplot（y1t,T，yt）； subplot(3，1，3);

stem(f,abs(fyt)/max(abs（fyt）)，’。’); grid;

title（’（c) s（t）的频谱’）; axis（[0，Fs/5,0,1。2]）;

xlabel（'f/Hz');ylabel（’幅度'）； clc；clear

%子程序

%产生信号程序 function st=mstg N=800

FS=10000;T=1/FS;TP=N＊T； t=0：T:(N-1)＊T;K=0：N—1; f=K/TP;

fc1=FS/10；％第一路调幅信号的载波频率fc1=1000HZ

fm1=fc1/10；％第一路调幅信号的调制信号频率为fm1=100hz. fc2=FS/20；%第二路调幅信号的载波频率fc2=500HZ

fm2=fc2/10；％第二路调幅信号的调制信号频率为fm2=50hz。

9

fc3=FS/40;%第三路调幅信号的载波频率fc3=250HZ

fm3=fc3/10;％第三路调幅信号的调制信号频率为fm3=25hz。 xt1=cos(2*pi*fm1*t).*cos（2*pi＊fc1＊t)；

xt2=cos（2＊pi＊fm2＊t）.*cos（2＊pi＊fc2＊t); xt3=cos（2＊pi＊fm3＊t).*cos（2＊pi*fc3*t); st=xt1+xt2+xt3；

fxt=fft(st，N)；％计算信号st的频谱.%绘图 subplot（2，1，1） plot（t,st); grid on；

xlabel（'t/s’);ylabel（'s（t）’)；

axis（[0,TP/8，min（st），max(st)］）； title（'(a）s(t)的波形’) subplot（2,1，2)

stem(f，abs（fxt)/max(abs（fxt)），’。')； grid on；title('(b） s(t）的频谱'） axis（[0，FS/5，0，1.2］)； xlabel(’f/HZ');ylabel(’幅度’)；

%损耗输出波形

function myplot（B,A)

[H,W］=freqz(B,A,1000）； m=abs（H);

plot(W/pi,20＊log10(m/max（m)））; grid on；

xlabel（’\\omega/\\pi’）; ylabel(’幅度（dB）'); axis(［0，1,-80,5］);

title（’(a) s(t）损耗函数曲线’）;

％滤波器输出波形 function tplot(xn,T,yn) n=0:length(xn)—1; t=n＊T;

plot(t，xn)； grid on； xlabel（'t/s');

ylabel（’y（n）'）;

axis(［0，t（end），min(xn),1.2*max(xn）]） title('(b） s(t)的波形'）；

10