数字信号处理实验46资料.docx
《数字信号处理实验46资料.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验46资料.docx（28页珍藏版）》请在冰豆网上搜索。
数字信号处理实验46资料
实验4离散系统的变换域分析
一、实验目的
1、熟悉对离散系统的频率响应分析方法；
2、加深对零、极点分布的概念理解。
二、实验原理
离散系统的时域方程为
其变换域分析方法如下：
频域:
系统的频率响应为:
Z域:
系统的转移函数为:
分解因式:
，
其中
和
称为零、极点。
三、预习要求
1.在MATLAB中，熟悉函数tf2zp、zplane、freqz、residuez、zp2sos的使用，其中：
[z，p，K]=tf2zp（num，den）求得有理分式形式的系统转移函数的零、极点；zplane（z，p）绘制零、极点分布图；h=freqz（num,den,w）求系统的单位频率响应；[r，p，k]=residuez（num，den）完成部分分式展开计算；sos=zp2sos（z，p，K）完成将高阶系统分解为2阶系统的串联。
2.阅读扩展练习中的实例，学习频率分析法在MATLAB中的实现；
3.编程实现系统参数输入，绘出幅度频率响应和相位响应曲线和零、极点分布图。
四、实验内容
求系统
的零、极点和幅度频率响应和相位响应。
解析：
【代码】
num=[0.05280.07970.12950.12950.7970.0528];
den=[1-1.81072.4947-1.88010.9537-0.2336];
[z,p,k]=tf2zp（num,den）;
disp（'零点'）;disp（z）;
disp（'极点'）;disp（p）;
disp（'增益系数'）;disp（k）;
figure
（1）
zplane（num,den）
figure
（2）
freqz（num,den,128）
【图形】
【结果】
零点
-1.5870+1.4470i
-1.5870-1.4470i
0.8657+1.5779i
0.8657-1.5779i
-0.0669+0.0000i
极点
0.2788+0.8973i
0.2788-0.8973i
0.3811+0.6274i
0.3811-0.6274i
0.4910+0.0000i
增益系数
0.0528
五、扩展练习
例1:
求下列直接型系统函数的零、极点，并将它转换成二阶节形式
解:
用MATLAB计算程序如下：
num=[1-0.1-0.3-0.3-0.2];
den=[10.10.20.20.5];
[z,p,k]=tf2zp（num,den）;
m=abs（p）;
disp（'零点'）;disp（z）;
disp（'极点'）;disp（p）;
disp（'增益系数'）;disp（k）;
sos=zp2sos（z,p,k）;
disp（'二阶节'）;disp（real（sos））;
zplane（num,den）
输入到“num”和“den”的分别为分子和分母多项式的系数。
计算求得零、极点增益系数和二阶节的系数：
零点:
0.9615
-0.5730
-0.1443+0.5850i
-0.1443-0.5850i
极点:
0.5276+0.6997i
0.5276-0.6997i
-0.5776+0.5635i
-0.5776-0.5635i
增益系数:
1
二阶节:
1.0000-0.3885-0.55091.00001.155200.6511
1.00000.288500.363001.0000-1.05520.7679
系统函数的二阶节形式为：
极点图见图:
图系统函数的零、极点图
例2:
差分方程
所对应的系统的频率响应。
解:
差分方程所对应的系统函数为：
用MATLAB计算的程序如下：
k=256;
num=[0.8-0.440.360.02];
den=[10.7-0.45-0.6];
w=0:
pi/k:
pi;
h=freqz（num,den,w）;
subplot（2,2,1）;
plot（w/pi,real（h））;grid
title（'实部'）
xlabel（'\omega/\pi'）;ylabel（'幅度'）
subplot（2,2,2）;
plot（w/pi,imag（h））;grid
title（'虚部'）
xlabel（'\omega/\pi'）;ylabel（'Amplitude'）
subplot（2,2,3）;
plot（w/pi,abs（h））;grid
title（'幅度谱'）
xlabel（'\omega/\pi'）;ylabel（'幅值'）
subplot（2,2,4）;
plot（w/pi,angle（h））;grid
title（'相位谱'）
xlabel（'\omega/\pi'）;ylabel（'弧度'）
图
实验5有限冲激响应数字滤波器设计
一、实验目的：
1、加深对数字滤波器的常用指标理解。
2、学习数字滤波器的设计方法。
二、实验原理：
低通滤波器的常用指标：
（1）通带边缘频率
；
（2）阻带边缘频率
；
（3）通带起伏
；
（4）通带峰值起伏
，
（5）阻带起伏
，最小阻带衰减
。
三、预习要求
1.在MATLAB中，熟悉函数fir1、kaiserord、remezord、remez的使用；
B=fir1（n,Wn,'high','noscale'）设计滤波器；
[n,Wn,beta,ftype]=kaiserord（f,a,dev）估计滤波器阶数；
[n,fo,ao,w]=remezord（f,a,dev,fs）计算等波纹滤波器阶数n和加权函数w（ω）；
B=remez（n,f,a）进行等波纹滤波器的设计。
2.阅读扩展练习中的实例，学习FIR滤波器的设计方法及其在MATLAB中的实现；
3.给出FIR数字滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。
数字滤波器有IIR和FIR两种类型，它们的特点和设计方法不同。
四、实验内容：
利用MATLAB编程,分别用窗函数法和等波纹滤波器法设计两种FIR数字滤波器，指标要求如下：
通带边缘频率：
，通带峰值起伏：
。
阻带边缘频率：
，最小阻带衰减：
。
●用窗函数法设计：
【解析】
调用函数[n,wn,bta,ftype]=kaiserord（f，a，dev，fs）
参数：
f=[0.3 0.45 0.65 0.8]为对应数字频率
a=[0 1 0]为由f指定的各个频带上的幅值向量，一般只有0和1表示；和f长度关系为（2*a的长度）—2=（f的长度）
devs=[0.01 0.1087 0.01]用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差，长度与a相等，计算公式：
阻带衰减误差=
，通带波动衰减误差=
fs缺省值为2HZ
【代码】
[n,wn,bta,ftype]=kaiserord（[0.30.450.650.8],[010],[0.010.10870.01]）;
%用kaiserord函数估计出滤波器阶数n和beta参数
h1=fir1（n,wn,ftype,kaiser（n+1,bta）,'noscale'）;
[hh1,w1]=freqz（h1,1,256）;
figure
（1）
subplot（2,1,1）
plot（w1/pi,20*log10（abs（hh1）））
grid
xlabel（'归一化频率w'）;ylabel（'幅度/db'）
subplot（2,1,2）
plot（w1/pi,angle（hh1））
grid
xlabel（'归一化频率w'）;ylabel（'相位/rad'）;
【图形】
●用等波纹法设计：
【解析】
调用函数[n,fpts,mag,wt]=remezord（f，a，dev） f=[0.3 0.45 0.65 0.8] a=[0 1 0]
dev=[0.01 0.1087 0.01]
其含义同函数[n,wn,bta,ftype]=kaiserord（f，a，dev，fs）中的参数相同。
【代码】
[n,fpts,mag,wt]=remezord（[0.30.450.650.8],[010],[0.010.10870.01]）;
%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w，使remez函数设计出的滤波器满足f、a及dev指定的性能要求
h2=remez（n,fpts,mag,wt）;%设计出等波纹滤波器
[hh2,w2]=freqz（h2,1,256）;
figure
（2）
subplot（2,1,1）
plot（w2/pi,20*log10（abs（hh2）））
grid
xlabel（'归一化频率w'）;ylabel（'幅度/db'）
subplot（2,1,2）
plot（w2/pi,angle（hh2））
grid
xlabel（'归一化频率w'）;ylabel（'相位/rad'）;
【图形】
五、扩展练习
例1:
用凯塞窗设计一FIR低通滤波器，通带边界频率
，阻带边界频率
，阻带衰减
不小于50dB。
解:
首先由过渡带宽和阻带衰减
来决定凯塞窗的N和
上图给出了以上设计的频率特性，（a）为N=30直接截取的频率特性（b）为凯塞窗设计的频率特性。
凯塞窗设计对应的MATLAB程序为：
wn=kaiser（30,4.55）;
nn=[0:
1:
29];
alfa=（30-1）/2;
hd=sin（0.4*pi*（nn-alfa））./（pi*（nn-alfa））;
h=hd.*wn';
[h1,w1]=freqz（h,1）;
或者：
b=fir1（29,0.4,kaiser（30,4.55））;
[h1,w1]=freqz（b,1）;
plot（w1/pi,20*log10（abs（h1）））;
axis（[0,1,-80,10]）;
grid;
xlabel（'归一化频率/p'）;
ylabel（'幅度/dB'）;
还可以使用[n,Wn,beta,ftype]=kaiserord（f,a,dev）函数来估计滤波器阶数等，得到凯塞窗滤波器：
fcuts=[0.30.5];%归一化频率omega/pi
mags=[10];
devs=[0.0510^（-2.5）];
[n,Wn,beta,ftype]=kaiserord（fcuts,mags,devs）;%计算出凯塞窗N，beta的值
hh=fir1（n,Wn,ftype,kaiser（n+1,beta）,'noscale'）;
freqz（hh）;
实际中，一般调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W（ω），调用函数remez可进行等波纹滤波器的设计，直接求出滤波器系数。
函数remezord中的数组fedge为通带和阻带边界
例2：
利用雷米兹交替算法设计等波纹滤波器，设计一个线性相位低通FIR数字滤波器，其指标为：
通带边界频率fc=800Hz，阻带边界fr=1000Hz，通带波动
阻带最小衰减At=40dB，采样频率fs=4000Hz。
解:
在MATLAB中可以用remezord和remez两个函数设计，其结果如图2，MATLAB程序如下：
fedge=[8001000];
mval=[10];
dev=[0.05590.01];
fs=4000;
[N,fpts,mag,wt]=remezord（fedge,mval,dev,fs）;
b=remez（N,fpts,mag,wt）;
[h,w]=freqz（b,1,256）;
plot（w*2000/pi,20*log10（abs（h）））;
grid;
xlabel（'频率/Hz'）;
ylabel（'幅度/dB'）;
所得图像如下所示：
实验6无限冲激响应数字滤波器设计
一、实验目的
1、掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法;
2、熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。
二、实验原理
在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：
1）利用buttord和cheb1ord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；
2）[num,den]=butter（N,Wn）（巴特沃斯）和[num,den]=cheby1（N,Wn）,[num,den]=cheby2（N,Wn）（切比雪夫1型和2型）可以进行滤波器的设计；
3）lp2hp，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；
4）使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；
5）利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。
三、预习要求
1.在MATLAB中，熟悉函数butter、cheby1、cheby2的使用，其中：
[num,den]=butter（N,Wn）巴特沃斯滤波器设计；
[num,den]=cheby1（N,Wn）切比雪夫1型滤波器设计；
[num,den]=cheby2（N,Wn）切比雪夫2型滤波器设计。
2.阅读扩展练习中的实例，学习在MATLAB中进行数字滤波器的设计；
3.给出IIR数字滤波器参数和滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。
四、实验内容
利用MATLAB编程，用脉冲响应不变法和双线性变换法设计一个数字带通滤波器，指标要求如下：
通带边缘频率：
，
，通带峰值起伏：
；
阻带边缘频率：
，
，最小阻带衰减：
。
【解析】
●巴特沃兹原型
1 双线性变换法
【代码】
ws1=2*8000*tan（0.3*pi/2）;
ws2=2*8000*tan（0.8*pi/2）;
wp1=2*8000*tan（0.45*pi/2）;
wp2=2*8000*tan（0.65*pi/2）;
ws=[ws1ws2];
wp=[wp1wp2];
Rp=1;Rs=40;
[N,Wn]=buttord（wp,ws,Rp,Rs,'s'）;
[num,den]=butter（N,Wn,'s'）;
[B,A]=bilinear（num,den,8000）;
[h,w]=freqz（B,A）;
f=w/pi*4000;
subplot（2,1,1）;
plot（f,20*log10（abs（h）））;
axis（[0,4000,-60,10]）;
grid
xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;
subplot（2,1,2）
plot（f,angle（h））;
grid
xlabel（'频率/Hz'）;ylabel（'相位'）;
【图形】
2 脉冲响应不变法
【代码】
fs=8000;
ws1=0.3*pi*fs;ws2=0.8*pi*fs;
wp1=0.45*pi*fs;wp2=0.65*pi*fs;
ws=[ws1ws2];wp=[wp1wp2];
Rp=1;Rs=40;
[N,Wn]=buttord（wp,ws,Rp,Rs,'s'）;
[num,den]=butter（N,Wn,'s'）;[B,A]=impinvar（num,den,8000）;
[h,w]=freqz（B,A）;
f=w/pi*4000;
subplot（2,1,1）
plot（f,20*log10（abs（h）））;
axis（[0,4000,-80,10]）;
grid;
xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;
subplot（2,1,2）
plot（f,angle（h））
grid;
xlabel（'频率/Hz'）;ylabel（'相位'）;
【图形】
fs=2;
ws1=2*fs*tan（0.3*pi/2）;
ws2=2*fs*tan（0.8*pi/2）;
wp1=2*fs*tan（0.45*pi/2）;
wp2=2*fs*tan（0.65*pi/2）;
ws=[ws1ws2];wp=[wp1wp2];Rp=1;Rs=40;
[N,Wn]=buttord（wp,ws,Rp,Rs,'s'）;
[num,den]=butter（N,Wn,'s'）;
[B,A]=bilinear（num,den,fs）;
[h1,w]=freqz（B,A）;
ws11=0.3*pi*fs;ws22=0.8*pi*fs;
wp11=0.45*pi*fs;wp22=0.65*pi*fs;
ws0=[ws11ws22];wp0=[wp11wp22];
[N,Wn]=buttord（wp0,ws0,Rp,Rs,'s'）;
[num,den]=butter（N,Wn,'s'）;
[B,A]=impinvar（num,den,fs）;
[h2,w]=freqz（B,A）;
figure
（1）
plot（w/pi,20*log10（abs（h1））,'-.',w/pi,20*log10（abs（h2））,'-'）;
axis（[0,1,-80,0]）;
grid;
xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;
●切必雪夫原型
1 双线性变换法
【代码】
ws1=2*8000*tan（0.3*pi/2）;
ws2=2*8000*tan（0.8*pi/2）;
wp1=2*8000*tan（0.45*pi/2）;
wp2=2*8000*tan（0.65*pi/2）;
ws=[ws1ws2];wp=[wp1wp2];Rp=1;Rs=40;
[N,Wn]=cheb1ord（wp,ws,Rp,Rs,'s'）;
[num,den]=cheby1（N,1,Wn,'s'）;[B,A]=bilinear（num,den,8000）;
[h,w]=freqz（B,A）;
f=w/pi*4000;
subplot（2,1,1）
plot（f,20*log10（abs（h）））;
axis（[0,4000,-60,10]）;
grid;
xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;
subplot（2,1,2）
plot（f,angle（h））;
grid;
xlabel（'频率/Hz'）;ylabel（'相位'）;
【图形】
2 脉冲响应不变法
【代码】
fs=8000;
ws1=0.3*pi*fs;ws2=0.8*pi*fs;
wp1=0.45*pi*fs;wp2=0.65*pi*fs;
ws=[ws1ws2];wp=[wp1wp2];
Rp=1;Rs=40;
[N,Wn]=cheb1ord（wp,ws,Rp,Rs,'s'）;
[num,den]=cheby1（N,1,Wn,'s'）;
[B,A]=impinvar（num,den,8000）;
[h,w]=freqz（B,A）;
f=w/pi*4000;
subplot（2,1,1）;
plot（f,20*log10（abs（h）））;
axis（[0,4000,-90,10]）;
grid;
xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;
subplot（2,1,2）;
plot（f,angle（h））;
grid;
xlabel（'频率/Hz'）;ylabel（'相位'）;
【图形】
fs=8000;
ws1=2*fs*tan（0.3*pi/2）;
ws2=2*fs*tan（0.8*pi/2）;
wp1=2*fs*tan（0.45*pi/2）;
wp2=2*fs*tan（0.65*pi/2）;
ws=[ws1ws2];wp=[wp1wp2];
Rp=1;Rs=40;
[N,Wn]=cheb1ord（wp,ws,Rp,Rs,'s'）;
[num,den]=cheby1（N,1,Wn,'s'）;
[B,A]=bilinear（num,den,fs）;
[h1,w]=freqz（B,A）;
ws11=0.3*pi*fs;ws22=0.8*pi*fs;
wp11=0.45*pi*fs;wp22=0.65*pi*fs;
ws0=[ws11ws22];wp0=[wp11wp22];
[N,Wn]=cheb1ord（wp0,ws0,Rp,Rs,'s'）;
[num,den]=cheby1（N,1,Wn,'s'）;
[B,A]=impinvar（num,den,fs）;
[h2,w]=freqz（B,A）;
figure
（1）
plot（w/pi,20*log10（abs（h1））,'-.',w/pi,20*log10（abs（h2））,'-'）;
axis（[0,1,-80,10]）;
grid;
xlabel（'频率/Hz'）
ylabel（'幅度/dB'）;
【结论】
双线性变换法通过将数字频率的取值范围从0到
p对应到模拟频率
的范围0到
的范围，也就对应于模拟域中所有可能的频率值。
双线性变换法不会出现频率混叠，但非线性关系却导致数字滤波器的频率响应不能逼真地模仿模拟滤波器的频率响应。
脉冲响应不变法通过选择满足设计要求的模拟滤波器冲激响应h（t）的采样值的数字脉冲响应h[n]得到的被采样的冲激响应将给出与原模拟滤波器非常相近的滤波器形状。
由于该方法不可避免的要发生频率混叠现象，所以只适合设计低通和带通滤波器。
从实验结果可以看出：
双线性变换法所设计的巴特沃兹滤波器最符合设计指标，而用脉冲响应不变法设计的滤波器（无论是巴特沃兹还是切比雪夫）有一定的误差，主要是由于混叠所引起。
五、扩展练习
例1：
设采样周期T=250μs（采样频率fs=4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc=1kHz。
[B,A]=butter（3,2*pi*1000,'s'）;
[num1,den1]=impinvar（B,A,4000）;
[h1,w]=freqz（num1,den1）;
[B,A]=butter（3,2/0.00025,'s'）;
[num2,den2]=bilinear（B,A,4000）;
[h2,w]=freqz（num2,den2）;
f=w/pi*2000;
plot（f,abs（h1）,'-.',f,abs（h2）,'-'）;
grid;
xlabel（'频率/Hz'）
ylabel（'幅值/dB'）
程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给
升级会员 