数字信号实验指导书带源程序.docx
《数字信号实验指导书带源程序.docx》由会员分享，可在线阅读，更多相关《数字信号实验指导书带源程序.docx（26页珍藏版）》请在冰点文库上搜索。
数字信号实验指导书带源程序
实验一离散时间系统与MATLAB
一.实验目的
1.进一步加深对离散时间系统的理解。
2.学习在MATLAB中怎样表示离散时间信号。
3.熟悉离散时间信号的作图。
二.实验步骤
1.复习离散时间系统的有关内容。
2.复习MATLAB的基本语法。
3.按实验内容熟悉stem。
4.编写程序。
5.输出结果，总结结论，按要求写出实验报告。
三.实验内容
1．掌握stem函数
STEM（Y）plotsthedatasequenceYasstemsfromthexaxisterminatedwithcirclesforthedatavalue.
STEM（X,Y）plotsthedatasequenceYatthevaluesspecifiedinX.
例：
t=[0:
0.1:
2];x=cos（pi*t+0.6）;stem（t,x）;
xn=[4,2,2,3,6,7];stem（xn）;
思考：
STEM（Y）与STEM（X,Y）有什么不同？
STEM与PLOT函数有什么不同？
2．掌握subplot函数
H=SUBPLOT（m,n,p）,orSUBPLOT（mnp）,breakstheFigurewindowintoanm-by-nmatrixofsmallaxes,selectsthep-thaxesforthecurrentplot,andreturnstheaxishandle.TheaxesarecountedalongthetoprowoftheFigurewindow,thenthesecondrow,etc.
例：
n1=0:
3;x1=[1,1,1,1];subplot（221）;stem（n1,x1）;title（'x1序列'）;
n2=0:
7;x2=[1,2,3,4,4,3,2,1];subplot（222）;stem（n2,x2）;title（'x2序列'）;
n3=0:
7;x3=[4,3,2,1,1,2,3,4];subplot（223）;stem（n3,x3）;title（'x3序列'）;
n4=0:
7;x41=cos（（pi/4）*n4）;subplot（224）;stem（n4,x41）;title（'x4序列'）;
思考：
subplot是怎样分配各个作图分区的顺序号的？
3．信号的运算
，
，请作出
，
的图形。
思考：
假如
与
长度不同，序列的求和和乘积运算能否执行，结果怎样？
程序：
x1=[1,0.7,0.4,0.1,0]
x2=[0.1,0.3,0.5,0.7,0.9]
subplot（121）;stem（x1+x2）;title（'x1+x2'）;
subplot（122）;stem（x1.*x2）;title（'x1*x2'）;
4．掌握freqs与freqz函数
freqs用于s域的频率响应的计算，freqz用于z域的频率响应的计算。
请参看两者的help文件。
例：
建立2个函数文件，分别输入为系统函数的分子，分母，和角频率的最大值，输出为系统频率响应的db，幅度，相位，和角频率。
参看下面的程序。
function[db,mag,pha,Omega]=freqs_m（b,a,Omega_Max）
%s域频率响应的计算
Omega=[0:
1:
500]*Omega_Max/500;
H=freqs（b,a,Omega）;
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
function[db,mag,pha,w]=freqz_m（b,a,w_Max）;
%z域频率响应的计算
w=[0:
1:
499]*w_Max/500;
H=freqz（b,a,w）;
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
（1）若
，请画出该系统幅频和相频特性。
程序：
b=1
a=[156]
Omega_Max=2*pi
[db,mag,pha,Omega]=freqs_m（b,a,Omega_Max）
subplot（221）;plot（Omega,mag）;title（'AF的幅度响应'）;
subplot（222）;plot（Omega,db）;title（'AF的幅度响应db'）;
subplot（223）;plot（Omega,pha）;title（'AF的相位响应'）;
（2）若
，请画出该的幅频和相频特性。
程序：
b=[11]
a=[10.80.64]
w_Max=2*pi
[db,mag,pha,w]=freqz_m（b,a,w_Max）;
subplot（221）;plot（w,mag）;title（'AF的幅度响应'）;
subplot（222）;plot（w,db）;title（'AF的幅度响应db'）;
subplot（223）;plot（w,pha）;title（'AF的相位响应'）;
思考：
freqs，freqz的输入的各个参数与系统函数怎么联系起来，各种调用形式之间有什么不同点？
观察上述系统的幅频和相频曲线，能得出什么结论？
四.实验报告要求
1.简述实验目的及原理。
2.概括各函数的常用调用方式，记录程序运行的各种结果。
3.回答思考题。
实验二用FFT作谱分析
一.实验目的
1.进一步加深DFT算法原理和基本性质的理解（FFT是DFT的一种快速算法，FFT的运算结果必然满足DFT的基本性质）。
2.熟悉FFT算法原理和FFT函数的应用。
3.学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因。
二.实验步骤
1.复习DFT的定义、性质和用DFT作谱分析的有关内容。
2.复习FFT算法原理与编程思想。
熟悉FFT算法的MATLAB实现。
注：
MATLAB提供fft函数来计算
的DFT，fft函数是用机器语言，而不是以MATLAB指令写成的，因此执行速度快。
格式
（1）：
y=fft（x）计算x的FFT变换y。
当x为矩阵，计算x中每一列信号的离散傅氏变换。
当x的长度为2的幂时，采用基2算法，否则采用分裂基算法。
格式
（2）：
y=fft（x,n）计算x的n点FFT，当x长度大于n时，截断x，否则补零。
Plot线性绘图函数。
stem：
绘制离散序列图。
subplot：
多坐标设置与定位当前坐标系。
figure：
创建新的图形窗口（用于输出图形的窗口）。
3.产生下列信号并进行谱分析。
（1）
（2）
（3）
（4）
（5）
（6）
4.编写程序
5.输出结果，总结结论，按要求写出实验报告。
三.实验内容
1.对上述6个信号，逐个进行谱分析
对于
，
，
，
，
：
。
对于
：
在一个周期内取
点抽样，做频谱分析。
（参考）每个信号谱分析的流程设计为：
程序：
%fft
figure
（1）;
n1=0:
3;x1=[1,1,1,1];subplot（221）;stem（n1,x1）;title（'x1序列'）;
k1=0:
7;y11=fft（x1,8）;magy11=[abs（y11）];subplot（222）;stem（k1,magy11）;title（'x1的8点FFT'）;
k2=0:
15;y12=fft（x1,16）;magy12=[abs（y12）];subplot（224）;stem（k2,magy12）;title（'x1的16点FFT'）;
%x2
figure
（2）;
n2=0:
7;x2=[1,2,3,4,4,3,2,1];subplot（221）;stem（n2,x2）;title（'x2序列'）;
k1=0:
7;y21=fft（x2,8）;magy21=abs（y21）;subplot（222）;stem（k1,magy21）;title（'x2的8点FFT'）;
k2=0:
15;y22=fft（x2,16）;magy22=abs（y22）;subplot（224）;stem（k2,magy22）;title（'x2的16点FFT'）;
%x3
figure（3）;
n3=0:
7;x3=[4,3,2,1,1,2,3,4];subplot（221）;stem（n3,x3）;title（'x3序列'）;
k1=0:
7;y31=fft（x3,8）;magy31=abs（y31）;subplot（222）;stem（k1,magy31）;title（'x3的8点FFT'）;
k2=0:
15;y32=fft（x3,16）;magy32=abs（y32）;subplot（224）;stem（k2,magy32）;title（'x3的16点FFT'）;
%x4
figure（4）;
n41=0:
7;x41=cos（（pi/4）*n41）;subplot（221）;stem（n41,x41）;title（'x4的8点序列'）;
k1=0:
7;y41=fft（x41,8）;magy41=abs（y41）;subplot（222）;stem（k1,magy41）;title（'x4的8点FFT'）;
n42=0:
15;x42=cos（（pi/4）*n42）;subplot（223）;stem（n42,x42）;title（'x4的16点序列'）;
k2=0:
15;y42=fft（x42,16）;magy42=abs（y42）;subplot（224）;stem（k2,magy42）;title（'x4的16点FFT'）;
%x5
figure（5）;
n51=0:
7;x51=sin（（pi/8）*n51）;subplot（221）;stem（n51,x51）;title（'x5的8点序列'）;
k1=0:
7;y51=fft（x51,8）;magy51=abs（y51）;subplot（222）;stem（k1,magy51）;title（'x5的8点FFT'）;
n52=0:
15;x52=sin（（pi/8）*n52）;subplot（223）;stem（n52,x52）;title（'x5的16点序列'）;
k2=0:
15;y52=fft（x52,16）;magy52=abs（y52）;subplot（224）;stem（k2,magy52）;title（'x5的16点FFT'）;
%x6
figure（6）;
f1=4;f2=8;f3=10;t=0:
0.01:
1;
x6=cos（f1*2*pi*t）+cos（f2*2*pi*t）+cos（f3*2*pi*t）;subplot（221）;plot（t,x6）;title（'x6模拟信号'）;
%x618
T=0.5/8;
t=0:
T:
0.5;
x6=cos（f1*2*pi*t）+cos（f2*2*pi*t）+cos（f3*2*pi*t）;
k61=0:
7;
y61=fft（x6,8）;magy61=abs（y61）;subplot（222）;stem（k61,magy61）;title（'x6的8点FFT'）;
%x6216
T=0.5/16;
t=0:
T:
0.5;
x6=cos（f1*2*pi*t）+cos（f2*2*pi*t）+cos（f3*2*pi*t）;
k62=0:
15;
y62=fft（x6,16）;magy62=abs（y62）;subplot（224）;stem（k62,magy62）;title（'x6的16点FFT'）;
2.令
，求出其
的FFT
。
当已求出
，根据DFT的对称性，求出
和
，与起第1步求出的
、
比较。
程序：
%x7
figure（7）;
n7=0:
15;x7=x42+x52;subplot（221）;stem（n7,x7）;title（'x7的16点序列'）;
k7=n7;y7=fft（x7,16）;magy7=abs（y7）;subplot（222）;stem（k7,magy7）;title（'x7的16点FFT'）;
X4=real（y7）;X5=imag（y7）*i;subplot（223）;stem（k7,abs（X4））;title（'x4的16点序列'）;
subplot（224）;stem（k7,abs（X5））;title（'x5的16点序列'）;
3.令
，求出其
的FFT
。
当已求出
，根据DFT的对称性，求出
和
，与起第1步求出的
、
比较。
程序：
%x8
figure（8）;
n8=0:
15;x8=x42+x52*j;subplot（221）;stem（n8,abs（x8））;title（'x8的16点序列'）;
k8=n8;y8=fft（x8,16）;magy8=abs（y8）;subplot（222）;stem（k8,magy8）;title（'x8的16点FFT'）;
tempk=mod（（16-k8）,16）;
tempy8=y8（tempk+1）;
convy8=real（tempy8）-imag（tempy8）*i;
X4=0.5*（y8+convy8）;subplot（223）;stem（k8,abs（X4））;title（'x4的16点序列'）;
X5=1/（2*j）*（y8-convy8）;subplot（224）;stem（k8,abs（X5））;title（'x5的16点序列'）;
四.思考题
1.在
时，
和
的幅频特性会相同吗？
为什么？
呢？
2.如果周期信号的周期预先不知道，如何使用FFT进行分析。
五.实验报告要求
1.简述实验的原理和目的。
2.输出实验中所用的程序以及所得的时域频域性曲线，与理论结果比较，分析误差产生的原因，以及FFT做谱分析时有关参数的选择方法。
3.总结实验所得的主要结论。
4.简要回答思考题。
实验三用双线性变换法设计IIR数字滤波器
一.实验目的
1.熟悉用双线性变换法设计IIR数字滤波器的原理与方法。
2.掌握数字滤波器的计算机仿真方法。
3.通过对实际心电图信号的滤波作用，获得数字滤波器的感性知识。
二.实验内容
人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。
下面给出一实际心电图信号采样序列样本
，其中存在高频干扰。
在实验中，以
作为输入序列，滤除其中的干扰成分。
要求：
程序设计应包含以下几个部分：
（1）绘制原始数据图形；
（2）设计巴特沃思低通滤波器并绘制其幅频相应曲线；（3）用设计的滤波器对原始数据进行滤波；（4）绘制滤波后的数据图。
三.实验步骤
1.用双线性变换法设计一个巴特沃思低通IIR数字滤波器。
设计指标参数为：
在通带内频率低于
时，最大衰减小于1dB，在阻带内
频率区间上，最小衰减大于15dB。
以
为间隔，输出模拟滤波器原型和数字滤波器在频率区间
上的幅频响应特性曲线。
设
。
2.用所设计的数字滤波器对实际心电图信号采样序列进行仿真滤波处理，输出比较滤波前后的心电图信号波形图，观察总结滤波作用。
3.（参考）可用到的matlab命令：
（1）文件I/O：
fopen（打开文件）fscanf（读入文件）fclose（关闭文件）
（2）信号处理内部函数：
buttap（模拟巴特沃思滤波器）poly（多项式系数）
freqs（模拟滤波器频响）freqz（数字滤波器频响）
filter（滤波）
（3）作图函数：
figure（新开图形窗口）subplot（图形定位）title（图形标题）
plot（线性图形）stem（序列图形）
函数的具体应用请参见help文件。
4.（参考）程序流程图：
wp=0.2*pi;wr=0.3*pi;Ap=1;Ar=15;%确定数字滤波器指标参数
T=1;Omegap=（2/T）*tan（wp/2）;Omegar=（2/T）*tan（wr/2）;%模拟滤波器指标
[cs,ds]=afd_butt（Omegap,Omegar,Ap,Ar）;%得到H（s）,其中cs为分子多项式系数,ds为分母多项式系数
[db,mag,pha,Omega]=freqs_m（cs,ds,pi）;%模拟滤波器的频率响应
figure;
subplot（231）;plot（Omega,mag）;title（'AF的幅度响应'）;%模拟滤波器响应
subplot（232）;plot（Omega,db）;title（'AF的幅度响应db'）;
subplot（233）;plot（Omega,pha）;title（'AF的相位响应'）;
[b,a]=bilinear（cs,ds,1/T）;%双线性变换法,从模拟滤波器到数字滤波器
[db,mag,pha,w]=freqz_m（b,a,pi）;%数字滤波器的频率响应
subplot（234）;plot（w/T,mag）;title（'DF的幅度响应'）;%数字滤波器的频率响应
subplot（235）;plot（w/T,db）;title（'DF的幅度响应db'）;
subplot（236）;plot（w/T,pha）;title（'DF的相位响应'）;
%开始滤波
fid=fopen（'text1.txt','r'）;
iffid==0
error（'cannotopentheinputfile!
'）
end
[inputx,N]=fscanf（fid,'%d'）;
figure;
n=0:
1:
（N-1）;
subplot（211）;stem（n,inputx）;title（'未滤波前的心电图'）;
y1=filter（b,a,inputx）;%y1=myfilter（inputx）;
subplot（212）;stem（n,y1）;title（'滤波后的心电图'）;
fclose（fid）;
function[db,mag,pha,Omega]=freqs_m（b,a,Omega_Max）
%s域频率响应的计算
Omega=[0:
1:
500]*Omega_Max/500;
H=freqs（b,a,Omega）;
mag=abs（H）;%幅度响应
db=20*log10（（mag+eps）/max（mag））;%幅度响应的db值
pha=angle（H）;%相位响应
function[db,mag,pha,w]=freqz_m（b,a,w_Max）;
%z域频率响应的计算
w=[0:
1:
499]*w_Max/500;
H=freqz（b,a,w）;
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
function[b,a]=afd_butt（Omegap,Omegar,Ap,Ar）;%求H（s）
%b分子多项式系数
%a分母多项式系数
%Omegap通带截止频率（rad/s）Ap通带衰减（dB）
%Omegar阻带截止频率（rad/s）Ar阻带衰减（dB）
krp=（（10^（0.1*Ar）-1）/（10^（0.1*Ap）-1））^0.5;
lemdarp=Omegar/Omegap;
N=ceil（log10（krp）/log10（lemdarp））;%计算滤波器阶数
Omegac=Omegar/（（10^（0.1*Ar）-1）^（1/2/N））;%3db截止频率通带满足指标，阻带指标有剩余
[b,a]=butter（N,Omegac,'s'）;
四.思考题
用双线性变换法设计数字滤波器过程中，变换公式
中
的取值，对设计结果有无影响？
为何？
五.实验报告要求
1.简述实验目的及原理。
2.由输出的
曲线及设计过程简述双线性变换的特点。
3.对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程和作用。
4.简要回答思考题。
实验四用窗函数法设计FIR数字滤波器
一.实验目的
1.掌握用窗函数法设计FIR数字滤波器的原理和方法
2.熟悉线性相位FIR数字滤波器特性
3.了解各种窗函数对滤波特性的影响
二.实验原理与方法
设所要求的理想数字滤波器的频率响应为
，
是与其对应的单位脉冲响应，因此：
由于
是矩形频率特性，故
是无限长的非因果序列。
而我们所要设计的是FIR数字滤波器，其单位脉冲响应
必然是有限长的，所以要用有限的
来逼近无限长的
，最有效的方法是截断
，用有限长的窗函数
来截取
，表示为：
这种方法即为窗口设计法。
三.实验内容和步骤
1.实验内容
（1）用海明窗设计一线性相位低通FIR数字滤波器，截止频率
。
窗口长度N=15、33。
要求在两种窗口长度情况下，分别求出
，输出相应的幅频特性和相频特性曲线，观察3dB和20dB带宽。
总结窗口长度N对滤波特性的影响。
%不同N对滤波器的影响
figure;
N=15;%窗口长度
wc=0.25*pi;%截止频率
n=0:
1:
N-1;
hd=ideal_lp（wc,N）;%求理想滤波器单位脉冲响应
subplot（221）;stem（n,hd）;title（'理想单位脉冲响应'）;
w_ham=（hamming（N））';
h=hd.*w_ham;%加窗
subplot（223）;stem（n,h）;title（'实际单位脉冲响应'）;
[H,w]=freqz（h）;%计算数字滤波器频响函数
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
subplot（222）;plot（w/pi,mag）;title（'实际的幅度响应'）;
subplot（224）;plot（w/pi,pha）;title（'实际的相位响应'）;
figure;
plot（w/pi,db）;title（'实际的幅度db响应'）;
db=20*log10（（mag+eps）/max（mag））;
figure;
N=33;
wc=0.25*pi;
n=0:
1:
N-1;
hd=ideal_lp（wc,N）;
subplot（221）;stem（n,hd）;title（'理想单位脉冲响应'）;
w_ham=（hamming（N））';
h=hd.*w_ham;
subplot（223）;stem（n,h）;title（'实际单位脉冲响应'）;
[H,w]=freqz（h）
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
subplot（222）;plot（w/pi,mag）;title（'实际的幅度响应'）;
subplot（224）;plot（w/pi,pha）;title（'实际的相位响应'）;
figure;
升级会员 