DSPB实验指导书.docx
《DSPB实验指导书.docx》由会员分享，可在线阅读，更多相关《DSPB实验指导书.docx（18页珍藏版）》请在冰点文库上搜索。
DSPB实验指导书
实验1Matlab的基本应用
一、实验目的
1、掌握MATLAB软件的基本使用方法。
2、掌握变量和文件名的命名规则以及M文件编辑器的基本使用方法。
3、学会使用MATLAB的帮助系统。
4、掌握数组的构造方法和数组、矩阵的四则运算。
5、掌握绘图常用函数的使用。
6、掌握简单的脚本文件和函数文件的编写及调试方法。
二、实验内容与要求
1.已知
，
，用MATLAB分别执行下列语句。
并在实验报告中记录语句和结果。
a.
b.
c.
d.
注意：
变量名的命名规则。
2.在M文件编辑器中编写程序画出函数
的图形，要求标注相关信息（如标题、横坐标、纵坐标等）以增加图形的可读性。
3.构造数组并回答以下指令运行后的结果是什么。
在实验报告中列出方法和结果
（1）求矩阵A的大小，（提示：
函数size）。
并回答它与length（A）的区别是什么？
（2）A（:
3）A（3,:
）A（2，3）A（2：
end,:
）A（1,[13]）
（3）列出数组A中所包含的数值0.6的下标。
提示：
用函数find。
（4）A（1：
2,[14]）=[2021;2223]
4.构造以下数组并回答以下表达式的结果是什么。
在实验报告中列出方法和结果
（1）r=a+b;
（2）r=a*d;（3）r=a.*d;
（4）r=a*c;（5）r=a.*c;（6）r=a\b;
（7）r=a.\b;（8）r=a.^b;（9）r=a.*[cc]
5.运行下列语句，回答以下命令运行的结果。
round（3.4）
ceil（3.4）
floor（3.4）
fix（3.4）
round（3.5）
ceil（3.5）
floor（3.5）
fix（3.5）
round（-3.4）
ceil（-3.4）
floor（-3.4）
fix（-3.4）
round（-3.5）
ceil（-3.5）
floor（-3.5）
fix（-3.5）
Mod（12,5）
Rem（12,5）
Mod（-12,5）
Rem（-12,5）
回答这几个函数的作用，并比较它们之间的区别。
6.编写脚本文件。
输入一组数据（数据可为正、负或零），根据以下公式计算其平均值和标准方差。
（注：
采用for循环）
注意：
学习程序的调试方法。
7.（选做）编写函数文件，通过输入不同的x,y值，根据以下公式计算f（x,y）的值并显示。
注意：
函数文件的编写格式和调用方法。
实验2序列的时域表示与运算
一、实验目的
1.掌握离散时间信号的表示；
2.掌握离散时间信号的基本运算（加、减、乘、反折、移位）的规则；
3.能用MATLAB进行简单的编程；
4.学习MATLAB函数的调用，实现序列的显示和运算。
二、实验内容与要求
请在实验报告中记录以下程序和结果
1．用MATLAB产生并画出下列序列的样本。
1）
2）
3）
4）
，式中
是在[-1，1]之间均匀分布的随机序列;
2.学习函数的调用。
设
，求解并画出下面序列。
3.已知序列
，在一个图形窗中分别画出该序列的实部、虚部、幅值和相位图。
注：
学习abs、angle、real、imag、subplot、title函数的使用。
4.若线性时不变系统的单位样值响应为
，输入序列
，求系统的输出
，并画出其波形图。
（思考：
你可以用几种方法来实现？
）
三、实验所用部分函数如下
1.单位冲激序列（信号）生成函数impseq
[x,n]=impseq（n0,n1,n2）
2.阶跃序列（信号）生成函数stepseq
[x,n]=stepseq（n0,n1,n2）
3.序列（信号）相加函数sigadd
[y,n]=sigadd（x1,n1,x2,n2）
以上为MATLAB没有，需外加入的函数（将相应函数拷贝到自己当前目录下）
4.正（余）弦生成函数sin、cos
y=sin（x）,y=cos（x）（注意:
x以弧度为单位）
5.随机序列生成函数rand，用法如：
Y=rand（n）生成n×n阶的均匀分布随机阵；
Y=rand（m,n）生成m×n阶的随机阵；
rand返回在[0，1]区间上的一个随机数；
将上面的rand写成randn则可以生成均值为0、方差为1的正态分布的随机变量。
6．全1矩阵生成函数ones（m,n）：
生成m×n阶全1矩阵
7．全0矩阵生成函数zeros（m,n）：
生成m×n阶全0矩阵
8．离散序列绘图函数stem
stem（y）以1、2、3…为横坐标， y为纵坐标画杆形图；
stem（x,y）以x为横坐标， y为纵坐标画杆形图（x与y数据个数必须一致）；
stem（…,’fill’）选项’fill’指定杆顶为实心，若无此选项则默认空心。
9．线性坐标平面绘图函数plot
用法与stem 类似，具体用法可查看MATLAB帮助
以上为MATLAB内置函数（在此仅为同学复习MATLAB提供）
四、参考程序
1.典型序列函数
1）单位阶跃函数（thestepsequence）
function[x,n]=stepseq（n0,n1,n2）
%Generatesx（n）=u（n-n0）;n1<=n,n0<=n2
%------------------------------------------
%[x,n]=stepseq（n0,n1,n2）
%
if（（n0n2）|（n1>n2））
error（'argumentsmustsatisfyn1<=n0<=n2'）
end
n=[n1:
n2];
%x=[zeros（1,（n0-n1））,ones（1,（n2-n0+1））];
x=[（n-n0）>=0];
2）单位冲激函数（theimpulsesequence）
function[x,n]=impseq（n0,n1,n2）
%Generatesx（n）=delta（n-n0）;n1<=n,n0<=n2
%----------------------------------------------
%[x,n]=impseq（n0,n1,n2）
%
if（（n0n2）|（n1>n2））
error（'argumentsmustsatisfyn1<=n0<=n2'）
end
n=[n1:
n2];
%x=[zeros（1,（n0-n1））,1,zeros（1,（n2-n0））];
x=[（n-n0）==0];
2.基本运算函数
1）加法函数
function[y,n]=sigadd（x1,n1,x2,n2）
%implementsy（n）=x1（n）+x2（n）
%-----------------------------
%[y,n]=sigadd（x1,n1,x2,n2）
%y=sumsequenceovern,whichincludesn1andn2
%x1=firstsequenceovern1
%x2=secondsequenceovern2（n2canbedifferentfromn1）
%
n=min（min（n1）,min（n2））:
max（max（n1）,max（n2））;%durationofy（n）
y1=zeros（1,length（n））;y2=y1;%initialization
y1（find（（n>=min（n1））&（n<=max（n1））==1））=x1;%x1withdurationofy
y2（find（（n>=min（n2））&（n<=max（n2））==1））=x2;%x2withdurationofy
y=y1+y2;%sequenceaddition
2）反折函数
以纵坐标轴为中心，左右对折，调用matlab自带函数fliplr。
function[y,n]=sigfold（x,n）
%implementsy（n）=x（-n）
%-----------------------------
%[y,n]=sigfold（x,n）
y=fliplr（x）;
n=-fliplr（n）;
3）乘法函数
function[y,n]=sigmult（x1,n1,x2,n2）
%implementsy（n）=x1（n）*x2（n）
%-----------------------------
%[y,n]=sigmult（x1,n1,x2,n2）
%y=producesequenceovern,whichincludesn1andn2
%x1=firstsequenceovern1
%x2=secondsequenceovern2（n2canbedifferentfromn1）
%
n=min（min（n1）,min（n2））:
max（max（n1）,max（n2））;%durationofy（n）
y1=zeros（1,length（n））;y2=y1;%initialization
y1（find（（n>=min（n1））&（n<=max（n1））==1））=x1;%x1withdurationofy
y2（find（（n>=min（n2））&（n<=max（n2））==1））=x2;%x2withdurationofy
y=y1.*y2;%sequencemultiplication
4）移位函数
function[y,n]=sigshift（x,m,n0）
%implementsy（n）=x（n-n0）
%-----------------------------
%[y,n]=sigshift（x,m,n0）
n=m+n0;
y=x;
五、编程注意事项
1）变量和文件名的命名规则：
以英文字母开头，由字母、数字和下划线组成，并且不要与MATLAB内置函数名称相同。
2）数组运算是点运算（对应位置的数据之间的加减乘除），注意小黑点；
实验三序列的频域表示与分析
一、实验目的
1.考察抽样间隔对信号频谱的影响；
2.掌握用FFT做谱分析的方法。
二、实验内容与要求（在预习报告中写出对信号进行傅立叶分析的程序，并思考下面的问题）
1.参考程序Lab3_1.m是对正弦信号的抽样，先改变抽样间隔（哪个参数？
），取若不同的值作观察图形变化，你能得出什么样的结论？
2.用FFT对模拟信号进行傅里叶分析
以频率fs对以下信号抽样N点
xa（t）=cos（at）+cos（bt）+cos（ct）
相应的参数是
a=2*pi*6500,b=2*pi*7000,c=2*pi*9000
fs=32000，N=16
对这N点序列作N点DFT，观察其幅频特性，如果
X=fft（x）
w是频率坐标向量，你可以考虑用stem（w,abs（X））,plot（w,abs（X））,plot（w,abs（X）,'*'）来显示，然后确定用哪种显示方式。
注意安排信号的时域、频域的显示。
1）对N=16点的序列作16点DFT（FFT），观察其幅频特性；
2）对N=16点的序列作M=256点DFT（FFT），这意味着在x后补了M-N个0，再观察幅频特性；
3）对N=256点的序列作256点DFT（FFT），观察其幅频特性；
4）改变fs，令其分别为24000，19000，18000，17000，16000，你看到了什么，做怎样的解释？
1）在实验报告中回答：
离散信号如何通过对连续信号抽样来获得，并按照以上要求写出对xa（t）抽样后的结果（即xn的表达式）；
2）在实验报告中回答：
如何对信号进行频谱分析，fft函数的用法？
3）在实验报告中回答：
对原时域信号补零后，其幅频特性有何变化？
4）在实验报告中回答：
高密度谱和高分辨率谱之间有何区别？
5）在实验报告中回答：
改变时域信号的采样率fs对信号的幅频特性会产生什么影响？
为什么？
在实际中选取fs的依据是什么？
三、参考程序
%ProgramLab3_1
%IllustrationoftheSamplingProcessintheTime-Domain
clf;
t=0:
0.0005:
1;
f=13;
xa=cos（2*pi*f*t）;
subplot（2,1,1）
plot（t,xa）;grid
xlabel（'Time,msec'）;ylabel（'Amplitude'）;
title（'Continuous-timesignalx_{a}（t）'）;
axis（[01-1.21.2]）
subplot（2,1,2）;
%inthelabproject,shouldrunfor4differentsamplingperiods
T=1./26;%抽样间隔
n=0:
T:
1;
xs=cos（2*pi*f*n）;
k=0:
length（n）-1;
stem（k,xs）;grid;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'Discrete-timesignalx[n]'）;
axis（[0（length（n）-1）-1.21.2]）
%ProgramLab3_2
%determinethespectrumofdiscrete-timesignalx[n]=cos（2*pi*500t）,thesamplingfrequencydenotesbyfs
fs=1500;M=256;
T=1/fs;n=0:
T:
（M-1）*T;
xa=cos（2*pi*500*n）;%discrete-timesignal
k=0:
（M-1）;
Xa1=abs（fft（xa,M））;&magnitudespectrum
plot（k*fs/（M-1）,Xa1）;
实验四FIR数字滤波器的设计
一、实验目的
1）掌握MATLAB中设计FIR数字滤波器的常用函数，并编写简单程序。
2）熟悉用FDATool进行FIR滤波器设计与分析的方法。
二、实验原理
用MATLAB提供的函数，设计FIR数字滤波器
三、实验内容与要求
1.用Kaiser窗设计一个FIR数字带阻滤波器，对模拟信号
xa（t）=cos（at）+cos（bt）+cos（ct）
a=2*pi*6500，b=2*pi*7000,c=2*pi*9000
滤波，要求滤去7000Hz的频率成分。
系统采样率为fs=32000Hz
这同我们第二次实验。
但采样点数应该比较大，可以用N=4096。
滤波器的Rp=0.25dB，As=50dB，过渡带宽可以用模拟频率（例如200Hz）也可以用数字频率指定。
还可以改变As（比如30dB）观察滤波效果。
注意：
A）如何从滤波器指标要求导出其它设计参数，如确定窗口类型、滤波器的阶数等；
B）如何验证设计得到的滤波器是否满足设计指标的语句；
试着直接用我们讲的MATLAB函数fir1（）、fir2（）进行设计
2.利用FDATool设计一个如上要求的FIR滤波器，对比使用不同的方法设计的滤波器。
四、参考程序：
1、FIR1Demo1.m用fir1函数设计高通滤波器例，各种窗口都试了一下；
2、FIR1Demo2.m用fir1函数设计多通带滤波器例；
3、FIR2Demo3.m用fir2函数设计多通带滤波器例；
4、FIRLsRemesDemo.m用firls和remez函数设计例，fir2也在这里做了对照。
1和4项所列的文件注释详细些。
FIR1Demo1.m
%用窗函数设计线性相位FIR滤波器
%教科书讲到的6种窗都用一下
%本例未涉及如何确定滤波器的设计参数
figNo=1;%显示用窗口号，接连用了3个
n=48%滤波器阶数
Wn=0.35;%截止频率
beta=0.1102*（60-8.7）;%Kaiser窗参数，假设阻带要求60分贝衰减（p.253）
%生成窗口矩阵，各窗函数看教科书pp.243-53
%先创建一个空矩阵，然后再把各窗函数列向量加进去
Windows=[];
Windows=[Windows,rectwin（n+1）];
Windows=[Windows,bartlett（n+1）];
Windows=[Windows,blackman（n+1）];
Windows=[Windows,hann（n+1）];
Windows=[Windows,hamming（n+1）];
Windows=[Windows,kaiser（n+1,beta）];
%接下来用一个循环完成用各窗口函数设计的滤波器
%并得到相应的传输函数向量，按列放在矩阵Hs中
%在这个循环中，win依次取矩阵Windows的各列
Hs=[];
forwin=Windows
b=fir1（n,Wn,'high',win）;%把high参数去掉就是低通滤波器
[h,w]=freqz（b,1）;
Hs=[Hs,h];
end
%现在用3种尺度来显示,请观察对比各窗口
figure（figNo）
freqzplot（Hs,w,'linear'）
figure（figNo+1）
freqzplot（Hs,w,'squared'）
figure（figNo+2）
freqzplot（Hs,w）%用分贝数显示
figure（figNo）%把第一个窗口推到前头
FIR1Demo2.m
%用窗函数设计线性相位FIR带通/带阻滤波器
n=48;
Wn=[0.350.65];%通带或阻带
%Hamingwindowbydefault
b=fir1（n,Wn,'stop'）;%再省去第三个参数'stop'就是带通
win=rectwin（n+1）;%矩形窗，宽度是滤波器阶数+1
bb=fir1（n,Wn,'stop',win）;
[H,w]=freqz（b,1,512）;
[HH]=freqz（bb,1,512）;
TH=[H,HH];
freqzplot（TH,w,'squared'）
FIR2Demo3.m
%用窗函数设计线性相位FIR滤波器
%设计2个带通/带阻滤波器演示
n=48
f=[0.20.40.60.8];%f的元素一定不能包含0和1！
bhm=fir1（n,f,'DC-1'）;%'DC-0'表示频率[0,0.2]是阻带，'DC-1'是通带
win=rectwin（n+1）;%上面用缺省海明窗，下面用的矩形窗
brec=fir1（n,f,'DC-1',win）;
%显示
[Hm,w]=freqz（bhm,1）;
[Hrec]=freqz（brec,1）;
H2=[Hm,Hrec];
freqzplot（H2,w,'linear'）
FIRLsRemesDemo.m
%最优化FIR线性相位滤波器设计例
%均方误差最小准则
%b=firls（n,f,m）
%最大误差最小化准则，雷米兹交换算法
%b=remez（n,f,m）
%参数：
%b:
返回的滤波器传输函数分子多项式的系数，也就是单位样本响应
%n:
滤波器的阶数
%f:
行向量，以0开始递增，各关键频率点，以1结束
%m:
行向量，维数与f相同，各关键频率点对应的响应幅度值
%n=20;%Filterorder
%f=[00.40.51];%Frequencybandedges
%m=[1100];%Desiredamplitudes
%下面的例子是用这两种方法设计的两个通带的滤波器
%再加上FIR2设计的作为对照
%useplot是控制显示方式的变量，参看下面显示部分
useplot=0;
%试着改动各参数观察一下，比如改变过渡带的宽度、阶数等
n=40;
m=[11001100];
f=[00.30.40.50.60.70.81];
b1=firls（n,f,m）;
b2=remez（n,f,m）;
%用FIR2设计
b3=fir2（n,f,m）;
%下面的代码显示比较它们的频率响应
[Hb1,w]=freqz（b1,1）;%分母多项式只有常数项1，没指定返回点数，缺省512点
[Hb2]=freqz（b2）;%这里把分母项省略了，返回的点数和上行一样，所以不再需要w
Hb3=freqz（b3）;
ifuseplot
%用一句plot画3条曲线，第2条用绿色，虚线，第3条红色；grid画出网格线
subplot（2,1,1）
plot（w/pi,abs（Hb1）,w/pi,abs（Hb2）,'g--',w/pi,abs（Hb3）,'r'）,grid
subplot（2,1,2）
%plot（w/pi,angle（Hb）,w/pi,angle（Hbb）,'r--'）,grid
plot（w/pi,unwrap（angle（Hb1））,w/pi,unwrap（angle（Hb2））,'g--',w/pi,unwrap（angle（Hb3））,'r'）
grid
else%也可以用下面的来显示
%Hb,Hbb是和w同维数的列向量（从Wokspace里可以看到）
%H是2列合成的矩阵
H=[Hb1,Hb2,Hb
升级会员 