IIR数字滤波器设计（数字信号处理）
一、实验目的
1.熟悉双线性变换法设计IIR数字滤波器的原理与方法。
2.掌握IIR数字滤波器的MATLAB实现方法，会调用ellipord()和ellip()
函数设计各种滤波器。
3.观察分析滤波器输入输出数据波形，理解数字滤波的概念。
二、实验原理及步骤
（一）实验原理-双线性变换法
数字滤波器是对数字信号实现滤波的线性时不变系统。数字滤波实质上是一种运算过程，实现对信号的运算处理。输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器本质上是一个完成特定运算的数字计算过程，也可以理解为是一台计算机。描述离散系统输出与输入关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。时域离散系统的频域特性:
其中
、
分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特性），
是数字滤波器的单位取样响应的频谱，又称为数字滤经过滤波后的频域响应。只要按照输入信号频谱的特点和处理信号的目的，适当选择
，使得滤波后
满足设计的要求，这就是数字滤波器的滤波原理。
数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应(IIR)数字滤波器和有限长冲激响应(FIR)数字滤波器。IIR 数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型来实现，其差分方程为：
系统函数为：
设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H(z)，使其频率响应H(z)满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。
基本设计过程如下：
1.将给定的数字滤波器指标转换成模拟滤波器的指标。
2.设计模拟滤波器。
3.将模拟滤波器转换成数字滤波器系统函数。
（二）实验步骤
1.在完成滤波器设计之后，采用filter()对输入信号进行滤波处理。调用信号产生函数mstg产生由抑制载波调制信号相加构成的复合信号st，该函数还会自动绘图显示其时域波形和幅频特性曲线，如下图1所示。
由图1中（a）和（b）可见，三路信号时域混叠无法在时域分离，但在频域是可以分离的，可以通过滤波的方法进行分离，即通过设计IIR滤波器，分离这三个不同频率的信号。
2.要求将三路信号进行分离，通过观察st信号的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。滤波器的通带最大衰减为0.1dB，阻带最小衰减为60dB。
3.调用ellipord()和ellip()分别设计这三个椭圆滤波器，并绘图显示其损耗函数曲线分别如图2，图3，图4。
图2 低通损耗函数曲线
4.调用filter()，用三个滤波器分别对信号产生函数mstg产生的信号进行滤波，分离出st中的不同载波频率信号，并绘图显示三个信号的时域波形，分别如图5，图6，图7。
三、实验结果分析
（一）函数mstg
阅读信号产生函数mstg，确定三路调幅信号的载波频率和调制信号频率。第1路调幅信号的载波频率fc1=1000Hz，调制信号频率fm1=100Hz。第2路调幅信号的载波频率fc2=500Hz，调制信号频率fm2=50Hz。第3路调幅信号的载波频率fc2=250Hz，调制信号频率fm2=25Hz。
- 采样点数对频谱图的影响
1.调幅信号
当N=800,1600,1800,2000时，调幅信号产生的频谱图如图8,9,10,11所示。
图8 N=800时s（t）的波形及频谱图
分析：因为信号s（t）是周期序列，频谱分析时要求观察时间为整数倍周期。s（t）的每个频率成分都是25Hz的整数倍。采样频率Fs=10kHz=25*400Hz，即在25Hz的正弦波的一个周期中采样400个点。所以，当N为400的整数倍时一定为s（t）的整数倍周期。因此，采样点数N=800,1600,2000时，对s（t）进行N点FFT可以得6根理想谱线，而当N=1800时，不是400的整数倍，则不能得到。
2.AM调幅信号
当N=800,1600,1800,2000时产生的频谱图如图12,13,14,15所示。
分析：因为信号s（t）时周期序列，频谱分析时要求观察时间为整数倍周期。因此，采样点数N=800,1600,2000时，对s（t）进行N点FFT可以得理想谱线，而当N=1800时，不是400的整数倍，则不能得到。当将该调幅信号修改为AM信号后，s（t）的频谱中有较大的频谱分量。如图所示。
- IIR滤波器
- 滤波器参数选取
由（一）可知，三路调幅信号的载波频率分别为250Hz，500Hz，1000Hz。带宽为50Hz，100Hz，200Hz。所以分离混合信号st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的指标参数选取如下：
对载波频率为250Hz的调幅信号，用低通滤波器分离，其通带截止频率fp=280Hz，通带最大衰减ap=0.1dB，阻带截止频率为fs=450Hz，阻带最小衰减as=60dB。对载波频率为500Hz的调幅信号，用低通滤波器分离，其通带截止频率fpl=440Hz、fph=560Hz，通带最大衰减ap=0.1dB，阻带截止频率为fsl=275Hz、fsh=900Hz阻带最小衰减as=60dB。对载波频率为1000Hz的调幅信号，用低通滤波器分离，其通带截止频率fp=890Hz，通带最大衰减ap=0.1dB，阻带截止频率为fs=600Hz，阻带最小衰减as=60dB。
IIR.m
代码语言：javascript复制function main
st=mstg; %调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号s
Fs=10000; T=1/Fs; %采样频率
fp=280; fs=450;wp=2*fp/Fs;ws=2*fs/Fs; %低通
name='y_1(t)';LHL_filter(wp,ws,name,st,T,'low');
fp_L=440; fp_R=560;fs_L=275; fs_R=900;%带通
wp=[2*fp_L/Fs,2*fp_R/Fs];ws=[2*fs_L/Fs,2*fs_R/Fs];
name='y_2(t)'; LHL_filter(wp,ws,name,st,T,'bandpass');
fp=890; fs=600;wp=2*fp/Fs;ws=2*fs/Fs; %高通
name='y_3(t)';LHL_filter(wp,ws,name,st,T,'high');
end
function LHL_filter(wp,ws,name,st,T,flag)
rp=0.1;rs=60; �指标（低通滤波器的通、阻带边界频）
[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp?
[B,A]=ellip(N,rp,rs,wp,flag); %调用ellip计算椭圆带通DF系统函数系数向量B和A
yt=filter(B,A,st); %滤波器软件实现
figure; freqz(B,A);%B,A分别为系统函数分子，分母多项式系数向量?
figure;tplot(st,yt,T,name);
end
function tplot(st,xn,T,name)
%时域序列连续曲线绘图函数?
%?xn:信号数据序列，%?T为采样间隔
N=length(xn);n=0:N-1; t=n*T; Tp=N*T; f=n/Tp;
subplot(311);plot(t,xn);xlabel('t/s');ylabel(name);
axis([0,t(end),min(xn),1.2*max(xn)]);
subplot(312);stem(f,abs(fftshift(fft(st,N))));
title('原 st(t)的频谱');xlabel('f/Hz');ylabel('幅度');
subplot(313);stem(f,abs(fftshift(fft(xn,N))));
title('滤波后xn(t)的频谱');xlabel('f/Hz');ylabel('幅度');
end
function st=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg 返回三路调幅信号相加形成的混合信号，长度N=1600
N=1600;Fs=10000;T=1/Fs;Tp=N*T; % N为信号st的长度。采样频率Fs=10kHz，Tp为采样时间
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率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(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/2,min(st),max(st)]);title('(a) s(t)的波形');
subplot(212);stem(f,abs(fxt)./max(abs(fxt)));
title('(频谱');xlabel('f/Hz');ylabel('幅度');
endmstg.m
代码语言：javascript复制function st=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg; 返回三路调幅信号相加形成的混合信号，长度N=1600
N=2000;Fs=10000;T=1/Fs;Tp=N*T; % N为信号st的长度。采样频率Fs=10kHz，Tp为采样时间
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率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);xlabel('t/s');ylabel('s(t)');
title(['N=',num2str(N),'s(t)的波形']);axis([0,Tp/8,min(st),max(st)]);
subplot(2,1,2);stem(f,abs(fxt)/max(abs(fxt)));
title(['N=',num2str(N),'s(t)的频谱']);xlabel('f/Hz');ylabel('幅度');
axis([0,Fs/5,0,1.2]);
endAM_mstg.m
代码语言：javascript复制function AM_mstg
N=2000; Fs=10000;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;fm1=fc1/10;%第1路调幅信号的载波频率fc1=1000Hz,调制信号频率fm1=100Hz
fc2=Fs/20;fm2=fc2/10;%第2路调幅信号的载波频率fc2=500Hz,调制信号频率fm2=50Hz
fc3=Fs/40;fm3=fc3/10;%第3路调幅信号的载波频率fc3=250Hz,调制信号频率fm3=25Hz
xt1=[1 cos(2*pi*fm1*t)].*cos(2*pi*fc1*t);
xt2=[1 cos(2*pi*fm2*t)].*cos(2*pi*fc2*t);
xt3=[1 cos(2*pi*fm3*t)].*cos(2*pi*fc3*t);
st=xt1 xt2 xt3;fxt=fft(st,N);
subplot(211);plot(t,st);xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title(['N=',num2str(N),' AM s(t)波形图'])
subplot(212);stem(f,abs(fxt)/max(abs(fxt)));title(['N= ',num2str(N),'频谱图'])
axis([0,Fs/5,0,1.2]);xlabel('f/Hz');ylabel('幅度');
end
