实验四IIR数字滤波器设计及软件实现实验报告.docx

实验四IIR数字滤波器设计及软件实现实验报告
一、实验目的熟悉用双线性变换法设计1IR数字滤波器的原理与方法；
(1) 学会调用MATLAB信号处理工具箱中滤波器设计函数(或滤波器设计分析工具 fdatool)设计各种IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。
(2) 掌握IIR数字滤波器的MATLAB实现方法。
(3)通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。
二、实验原理
设计11R数字滤波器一般采用间接法(脉冲响应不变法和双线性变换法),应用最广泛 的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器 的指标；②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函 数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六 章介绍的滤波器设计函数butter、chebyl、cheby2和ellip可以分别被调用来直接设计 巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函 数直接设计IIR数字滤波器。
本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给 定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n)。
三、实验内容及步骤
(1)调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st,该函 数还会自动绘图显示st的时域波形和幅频特性曲线，如图所示。由图可见，三路 信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离， 这就是本实验的目的。

(b) s(t)的频谱
•I•••4IIIII
| ------- ♦♦-♦- -- -4--'iiiii
翅；：；；：
los 0.514- — !1-e-IIIII
□200400600800 1000 1200 U00 1600 1800 2000f/Hz
图三路调幅信号st的时域波形和幅频特性曲线
要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分 第1页共12页
3. 具体步骤：
(1) 构造原始信号s(t)
(2) 画出s(t)的频谱
(3) 设计ellipse数字滤波器(IIR),包括低通，带通，带通，并显示幅频特 性
(4) 用得到的滤波器进行滤波，分离出三路信号，观察时域波形和幅频特性
(5) 用三路信号si, s2, s3尝试重新合成原始信号问题：为什么重新合成的信号和原信号不相等呢？谁能解释一下？谢谢 程序如下：
clearelf
% (1)构造原始信号Fs二10000; T二1/Fs; %先设定采样频率
t=0:T:0. 1;n=length(t);s=cos (2*pi*250*t). *cos(2*pi*25*t)+cos(2*pi*500*t). *cos(2*pi*50*t) +... cos (2*pi*1000*t)・ *cos(2*pi*100*t);
subplot (2, 1, 1), plot (t, s), axis ([0 0. 08 ~2 3])titleC原始信号s(t)')
xlabel (' t/s'), ylabel (' s (t)')% (2)画出s(t)的频谱
ft=fftshift(fft(s));i二fix(n/2) ;f=(-i:i)/n*Fs;%貌 似 这 是 公 式。
subplot(2, 1, 2), stem(f, abs(ft), 'Marker',' none'), xlim([0 1250]) title C s (t)的频谱')xlabel (J f/Hz，), ylabel C 幅度’)% ( 3 )设计 ellipse数字滤波器(IIR),并显示幅度特性%%3. la设计模拟低通滤 波 器 fp=320:fs=400:Ap=0. 1:As=60:
wp二2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp= (2*Fs)*tan(wp/2);Ws=(2*Fs)*tan(ws/2);
[N, Wc]= ellipord(Wp , Ws , Ap , As ,'s');[bLPs , aLPs]=ellip(N , Ap, As, Wc,' s');
[H, w]=freqs (bLPs, aLPs);db =20*logl0(abs(H)):
figure, subplot (2, 1, 1), plot (w/2/pi, db);axis ([0 1600 -80 5]),grid
title ('模拟低通滤波器的幅度特性');xxxxxx C f (Hz)');
yyyyyy (，dB')；%%3. lb将模拟低通滤波器转换为数字低通滤波器
[bLPz, aLPz] = bilinear (bLPs , aLPs , Fs);w 二 linspace (0, pi, 1000);
h = freqz(bLPz, aLPz , w);subplot (2,1, 2), plot (w*Fs/2/pi, 20*logl0(abs(h)));
axis ([0 1600 -80 5]),grid titlef数字低通滤波器的幅度特性');xxxxxx ('f (Hz)'); yyyyyyCdB')；
%%3. 2a设计模拟带通滤波器 fp=[430 570];fs=[330 670];Ap=0. l;As=60;wp二2*pi*fp/Fs;ws二2*pi*fs/Fs;
Wp= (2*Fs)*tan(wp/2);Ws=(2*Fs)*tan(ws/2);[N, Wc]= ellipord(Wp , Ws , Ap , As ,' s');
[bBPs , aBPs]=ellip (N , Ap, As, Wc,' s');[H, w] =freqs (bBPs, aBPs);
db =20*logl0(abs(H)): figure, subplot (2, 1, 1), plot (w/2/pi, db);axis ([0 1600 -80 5]),grid title ('模拟带通滤波器的幅度特性')；
xxxxxx (' f (Hz)');yyyyyy (' dB')
%%3. 2b将模拟带通滤波器转换为数字带通滤波器[bBPz,aBPz] = bilinear (bBPs , aBPs , Fs);
w = linspace (0, pi, 1000);h = freqz (bBPz, aBPz , w); subplot (2, 1, 2), plot(w*Fs/2/pi, 20*logl0(abs(h)));
axis([0 1600 -80 5]),grid titleC数字低通滤波器的幅度特性')； xxxxxx ('f (Hz)');yyyyyy ('dB')； %%3. 3a设计模拟高通滤波器
fp=800;fs=700;Ap=0. 1;As=60; wp二2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp=(2*Fs)*tan(wp/2);Ws=(2*Fs)*tan(ws/2);
[N, Wc]= ellipord (Wp , Ws , Ap , As ,'s');[bHPs , aHPs]=ellip (N, Ap, As, Wc ,' high',' s');
[H, w]=freqs (bHPs, aHPs);db 二20*logl0(abs(H)); figure, subplot (2, 1, 1), plot (w/2/pi, db);
axis([0 1600 -80 5]),grid titleC模拟高通滤波器的幅度特性')； xxxxxx (' f (Hz)');yyyyyy (，dB')
%%3. 3b将模拟高通滤波器转换为数字高通滤波器[bHPz, aHPz] = bilinear(bHPs , aHPs , Fs);
w = linspace(0, pi, 1000);h = freqz (bHPz, aHPz , w);
subplot (2, 1, 2), plot(w*Fs/2/pi, 20*logl0(abs(h))); axis([0 1600 -80 5]),grid title ('数字高通滤波器的幅度特性')；xxxxxx ('f (Hz)'); yyyyyy('dB')；
% (4)用得到的滤波器进行滤波，分离出三路信号，观察时域波形和幅频特 性％%4. 1滤波si = filter(bLPz, aLPz, s); s2 二 filter(bBPz, aBPz, s); s3 = filter (bHPz, aHPz, s);
%%4. 2时域波形figure, subplot (3, 1, 1), plot (t, si) xlabel (' t/s' ), ylabel (' si (t)') subplot (3, 1, 2), plot (t, s2) xlabel (' t/s'), ylabel (' s2 (t)') subplot (3, 1, 3), plot (t, s3) xlabel C t/s，), ylabel (' s3 (t)') %%4. 3幅频特性
ftl=fftshift(fft(si)):ft2=fftshift(fft(s2)) ;ft3=fftshift(fft(s3)); figure, subplot (3, 1, 1), stem(f, abs (ftl),' Marker5,' none'), xlim([0 1250]) subplot(3, 1, 2), stem(f, abs(ft2),' Marker5,'none'),xlim([0 1250]) subplot(3, 1, 3), stem(f, abs(ft3),' Marker5,' none，), xlim([0 1250]) %(5)用三路信号sl,s2, s3尝试重新合成原始信号ss=sl+s2+s3;%合成信号记作ss isequal (s, ss)%判断合成信号是否和原信号相等，结果不相等，why?
[本帖最后由cwjy于2010-5-318:19编辑]
离st中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器) 的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0. IdB,阻带最小衰减为60dB。
提示：抑制载波单频调幅信号的数学表示式为1
s(A) cos(2 ft) cos(2 ft) —[cos(2 (f f cos(2 (f f)/)]0c 2c oc 0
其中，cos(2々)称为载波，匚为载波频率，cos(2 Q)称为单频调制信号，f°为调制正 弦波信号频率，且满足。由上式可见，所谓抑制载波单频调幅信号，就是2个正弦信号相乘，它有2个频率成分：和频r f和差频f /,这2个频率成分关于载波频
c 0C 0
率 对称。所以，1路抑制载波单频调幅信号的频谱图是关于载波频率f对称的2根谱线， f CC其中没有载频成分，故取名为抑制载波单频调幅信号。容易看出，图10. 4.1中三路调幅信 号的载波频率分别为250Hz、500Hz、1000Hzo如果调制信号m(t)具有带限连续频谱，无直
流成分，则s(，) m(r)cos(2 ft)就是一般的抑制载波调幅信号。其频谱图是关于载波频率匚对称的2个边带(上下边带)，砖业课通信原理中称为双边带抑制载波 (DSB-SC) 调幅信号，简称双边带(DSB)信号。如果调制信号rn(t)有直流成分，则
s(f) m(0cos(2 ffl就是一般的双边带调幅信号。其频谱图是关于载波频率匚对称的2 个边带(上下边带),并包含载频成分。
(3) 编程序调用MATLAB滤波器设计函数ellipord和ellip分别设计这三个椭圆滤波 器，并绘图显示其幅频响应特性曲线。
(4) 调用滤波器实现函数filter,用三个滤波器分别对信号产生函数mstg产生的信 号st进行滤波，分离出st中的三路不同载波频率的调幅信号yjn)、y2(n)和y3(n),并 绘图显示yl(n)、y2(n)和y3(n)的时域波形，观察分离效果。
四、信号产生函数mstg清单
function st=mstg
%产生信号序列向量st,并显示st的时域波形和频谱 %st=mstg返回三路调幅信号相加形成的混合信号，长度N=1600
N=1600%N为信号st的长度。
Fs=10000;T=l/Fs;Tp=N*T; %采样频率 Fs=10kHz, Tp 为采样时间 t=O:T: (N-l)^T;k=O:N-l;f=k/Tp;
fcl=Fs/10;%第1路调幅信号的载波频率fcl=1000Hz,
fml=fcl/10;%第1路调幅信号的调制信号频率fml=100Hz
fc2二Fs/20;%第2路调幅信号的载波频率fc2=500IIz
fm2=fc2/10;%第2路调幅信号的调制信号频率fm2=50Hz
fc3=Fs/40;%第3路调幅信号的载波频率fc3=250Hz,
fm3=fc3/10;%第3路调幅信号的调制信号频率fm3=25Hz
xtl=cos(2*pi*fml*t)・ *cos(2*pi*fcl*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=xtl+xt2+xt3; fxt=fft (st,N);

st=xtl+xt2+xt3; fxt=fft (st,N);
%三路调幅信号相加 %计算信号st的频谱
%==以下为绘图部分，绘制st的时域波形和幅频特性曲线二 subplot (3, 1, 1)
plot (t, st) ;grid;xlabel (' t/s，) ;ylabel (' s (t)');axis([0, Tp/8, min (st), max (st)]) ; title (a) s(t)的波形') subplot (3, 1, 2)
stem(f, abs (fxt)/max(abs (fxt)),'・');grid; title(' (b) s(t)的频谱') axis ([0, Fs/5, 0, 1.2]); xlabel (' f/Hz') ;ylabel ('幅度')
五、实验程序框图如图io. 4. 2所示调用函数mstg产生st,自动绘图 显示st的时域波形和幅频特性曲线

调用ellipord和ellip分别设计三个椭圆
滤波器，并绘图显示其幅频响应特性曲线。
出三路不同载波频率的调幅信号y】（n）、y 2（n）和fn）

六、滤波器参数及实验程序清单
1、滤波器参数选取
观察图可知，三路调幅信号的载波 频率分别为250Hz、500Hz. lOOOHzo带宽
（也可以由信号产生 函数mstg清单看出）分别为50Hz、100Hz、200Hzo所以，分离混合 信号st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波 器）的指标参数选取如下：
对载波频率为250Hz的条幅信号，可以用低通滤波器分离，其指标为
带截止频率fp280Hz＞通带最大衰减p 0.1dBdB;
阻带截止频率f 450Hz＞阻带最小衰减s 60dBdB,
对载波频率为500Hz的条幅信号，可以用带通滤波器分离，其指标为
带截止频率。/440也，fpu 56C)hz,通带最大衰减p O.ldBdB;
阻带截止频率'275Hz)f 900Hz> hz,阻带最小衰减si
s 60dBdB,
对载波频率为1000Hz的条幅信号，可以用高通滤波器分离，其指标为
带截止频率。89°Hz,通带最大衰减p °」dBdB；
阻带截止频率55C)hz,阻带最小衰减、60dBdB,S
说明：(1)为了使滤波器阶数尽可能低，每个滤波器的边界频率选择原则是尽量使滤波 器过渡带宽尽可能宽。
(2) 与信号产生函数mstg相同，采样频率Fs=10kHzo为了滤波器阶数最低，选用椭圆滤波器。
按照图所示的程序框图编写的实验程序为exp4.ni。
2、实验程序清单
%实验4程序
% HR数字滤波器设计及软件实现
clear all;close all
Fs = 10000; T=l/Fs;%采样频率
%调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st
st=mstg;
%低通滤波器设计与实现
fp=280;fs=450;
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
ylt二filter (B, A, st) ;%滤波器软件实现
%低通滤波器设计与实现绘图部分
figure (2);subplot (3,1,1);
myplot(B, A); %调用绘图函数myplot绘制损耗函数曲线
yt- y_l(t)';
subplot (3, 1,2); tplot (ylt, T, yt) ; %调用绘图函数tplot绘制滤波器输出波形
%带通滤波器设计与实现
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, 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) ;%滤波器软件实现
%带通滤波器设计与实现绘图部分(省略)
%高通滤波器设计与实现
fp=890;fs=600;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0. l;rs=60;%DF指标(低通滤波器的通、阻带边界频)
[N, wp]=ellipord(wp, ws, rp, rs) ;%调用ellipord计算椭圆DF阶数N和通带截止频率wp
[B, A]=ellip(N, rp, rs, wp,' high，) ; %调用ellip计算椭圆带通DF系统函数系数向量B和A y3t=filter (B, A, st) ;%滤波器软件实现
%高低通滤波器设计与实现绘图部分(省略)
七、实验程序运行结果
实验4程序exp4.m运行结果如图104.2所示。由图可见，三个分离滤波器指标参数选 取正确，算耗函数曲线达到所给指标。分离出的三路信号yl(n), y2(n)和y3(n)的波形是抑 制载波的单频调幅波。
损耗函数曲线

(a)低通滤波器损耗函数及其分离出的调幅信号y (t)
i
损耗函数曲线
00.10.203040.506070.80.91

(b)带通滤波器损耗函数及其分离出的调幅信号y2(t)
O
00 010.Q20 030 040 050.060.07
t/§
（c）高通滤波器损耗函数及其分离出的调幅信号y3（t）
图104.实验4程序exp4.m运行结果
八、思考题简答
（1）请阅读信号产生函数mstg,确定三路调幅信号的载波频率和调制信号频率。（2） 信号产生函数nistg中采样点数N=800,对st进行N点FFT可以得到6根理想谱 线。如果取N二1000,可否得到6根理想谱线？为什么？ N二2000呢？请改变函数mstg中采 样点数N的值，观察频谱图验证您的判断是否正确。
（3）修改信号产生函数mstg,给每路调幅信号加入载波成分，产生调幅（AM）信号， 重复本实验，观察AM信号与抑制载波调幅信号的时域波形及其频谱的差别。
答：分析发现，st的每个频率成分都是25IIz的整数倍。采样频率Fs=10kHz=25X400IIZ, 即在25Hz的正弦波的1个周期中采样400点。所以，当N为400的整数倍时一定为st的 整数个周期。因此，采样点数N=800和N=2000时，对st进行N点FFT可以得到6根理想谱 线。如果取N=1000,不是400的整数倍，不能得到6根理想谱线
Mat lab的IIR滤波器描述清楚iir滤波器的设计建模过程，程序注释，以及对不同设计方案的说明比较等。
lirl :
低通巴特沃斯模拟滤波器设计。
通带截至频率3400 Hz ,通带最大衰减3dB阻带截至频率4000 Hz ,阻带最小衰减40dB
Iir2 :
模拟低通滤波器转换为数字低通滤波器,脉冲响应不变法和双线性变换法。
Iir3 :
切比雪夫二型低通数字滤波器设计通带边界频率0.2兀,通带最大衰减IdB
阻带截至频率0.47T ,阻带最小衰减80dBIir4 :
椭圆带通数字滤波器设计Iir5 :
高通和带通巴特沃思数字滤波器设计
双线性变换%低通巴特沃斯模拟滤波器设计
clear; close allfp=3400; fs=4000; Rp=3; As=40;
[N,fc]=buttord(fjp,fs,Rp,As」s)[B,A]=butter(N,fc,'s');
[hf,f] =freqs(B, A, 1024);plot(f,20*log 10(abs(hf)/abs(hf( 1))))
grid, xlabel('f/Hz'); ylabel('幅度(dB)') axis([0,4000,・40,5]);line([0,4000],[-3,-3]);
line([3400,3400],[-90,5])%用脉冲响应不变法和双线性变换法将模拟滤波器离散化 clear; close all
b=1000;a=[l,1000];w=[0:1000*2*pi];
[hf,w]=freqs(b9a,w);subplot(2,3,l); plot(w/2/pi,abs(hf)); grid;
xlabel('f(Hz)'); ylabel('幅度');title(，模拟滤波器频响特性)FsO=[ 1000,500];
for m=l :2Fs=FsO(m)
[d,c]=impinvar(b,a,Fs)[f,e]=bilinear(b,a,Fs)
wd=[0:512] *pi/512;hw 1 =freqz(d,c,wd);
h w2=freqz(f,e, wd);subplot(2,3,2); plot(wd/pi,abs(hwl)/abs(hwl(l))); grid on; hold on title(，脉冲响应不变法)
subplot(2,3,3); plot(wd/pi,abs(hw2)/abs(hw2( 1))); grid on; hold on title(，双线性变换法)end
%切比雪夫II型低通数字滤波器设计clear; close all
wp=0.2; ws=0.4; Rp=l; Rs=80;[N,wc]=cheb2ord(wp,ws,Rp,Rs)
[B,A]=cheby2(N,Rs,wc)freqz(B,A)
%直接设计带通数字椭圆滤波器clear; close all
Wp=[0.25,0.45]; Ws=[O.15055];Rp=0.1; Rs=60;
[N,wc]=ellipord(Wp,Ws,Rp,Rs)[b,a]=ellip(N,Rp,Rs,wc)
[hw,w]=freqz(b,a);subplot(2,l,l); plot(w/pi,20*logl0(abs(hw))); grid
axis([0,l,・80,5]); xlabel('w/ 兀')；ylabel('幅度(dB)')subplot(2,l,2); plot(w/pi,angle(hw)); grid
axis([0,l,-pi,pi]); xlabel(*w/ 兀')；ylabel('相位(rad)1)%用双线性变换法设计数字高通和带通滤波器 clear; close all
T=l; wch=pi/2;wlc=0.35*pi; wuc=0.65*pi;
B=1;A=[1,2.6131,3.4142,2.6131,1]; [h,w]=freqz(B,A,512);subplot(2,2,l); plot(w,20*logl0(abs(h))); grid
%axis([0,10,・90,0]); xlabel(w/ 兀');title(模拟低通幅度(dB)f) % 高通omegach=2*tan(wch/2)/T;[B hs, Ahs]=lp2hp(B, A,omegach); [Bhz,Ahz]=bilinear(Bhs,Ahs, 1/T);
[h,w]=freqz(Bhz, Ahz,512);subplot(2,2,3); plot(w/pi,20*logl0(abs(h))); grid
axis([0,1,-150,0]); xlabel('w/ n '); title('数字高通幅度(dB)') %带通omegalc=2*tan(wlc/2)/T; omegauc=2*tan(wuc/2)/T;
wo=sqrt(omegalc*omegauc); Bw=omegauc-omegalc; [Bbs,Abs]=lp2bp(B,A,wo,Bw);[Bbz,Abz]=bilinear(Bbs,Abs, 1/T);
[h,w]=freqz(Bbz,Abz,512);subplot(2,2,4); plot(w/pi,20*logl0(abs(h))); grid
axis([0,1,-150,0]); xlabel('w/ 兀')；title('数字带通幅度(dB)')数字滤波器设计与应用问题
1. 题目：数字滤波器的设计与应用 设计要求：利用Matlab软件，以复合信号分离为例，对“数字信号处理” 课程中的谱分析、数字滤波器设计和信号滤波这三个过程进行了仿真实现，给 出了仿真结果。