数字信号处理B实验指导书讲解.docx
《数字信号处理B实验指导书讲解.docx》由会员分享，可在线阅读，更多相关《数字信号处理B实验指导书讲解.docx（16页珍藏版）》请在冰点文库上搜索。
数字信号处理B实验指导书讲解
实验1矩阵、序列及系统的时域表示与运算
一、实验目的
掌握用MATLAB表示离散时间信号和系统，以及它们的运算和显示。
二、实验内容与要求
1．用MATLAB产生并画出下列序列的样本。
a．
b．
c．
d．
，式中
是在【-1，1】之间均匀分布的随机序列，你如何表征这个序列？
e．
为周期的，画出5个周期。
2若线性时不变系统的单位样值响应为
，输入序列
，求系统的输出
，并画出其波形图。
（思考：
你可以用几种方法来实现？
）
三、实验所用部分函数如下
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’指定杆顶为实心，若无此选项则默认空心。
stem（...,LineSpec）参数LineSpec指定杆形图的线形、数据标志符号及颜色，具体用法可查看MATLAB帮助
9．线性坐标平面绘图函数plot
用法与stem 类似，具体用法可查看MATLAB帮助
以上为MATLAB内置函数（在此仅为同学复习MATLAB提供）
四、在MATLAB程序中变量赋值注意问题
在MATLAB 中，对变量赋值时其维数可以按需要动态地改变，这样虽然方便程序设计但同时容易出错。
另外，频繁分配变量空间会大大降低程序的执行速度，因而应该尽量避免不必要的矩阵、向量维数的改变。
通常先用zeros（）函数给变量分配足够大小的空间，再对变量进行赋值。
例：
依次执行下面的语句
tic%开始计时
fori=1:
10000
c（i）=i;%每次都重新分配空间
end
toc%读取计时时间
tic%开始计时
d=zeros（1,10000）;%预先分配空间
fori=1:
10000
d（i）=i;%直接赋值，不必重新分配空间
end
toc%读取计时时间
运行结果如下：
elapsed_time=
1．1560
elapsed_time=
0．0470
从结果可以看出，第2种赋值方法所用的时间比第1种方法所用时间少得多（以上是在主频为2．66GHZ的机器上运行的结果）。
实验2离散信号、系统的频域表示
一、实验目的
1.考察抽样间隔对信号频谱的影响；
2.掌握用FFT做谱分析的方法。
二、实验内容与要求
1、用DFT/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）,'*'）来显示，然后确定用哪种显示方式。
接下来，对x作M=256点FFT，这意味着在x后补了M-N个0，再观察幅频特性；
现在令抽样点数N=256，再对这个抽样序列作N点FFT，观察幅频特性。
通过这些观察，你能作出什么判断？
试着改变fs，令其分别为24000，19000，18000，17000，16000，你看到了什么，做怎样的解释？
注意安排你的时域、频域的显示。
MATLAB程序如下
fs=[32000,24000,19000,18000,17000,16000];
m=length（fs）;
a=2*pi*6500;
b=2*pi*7000;
c=2*pi*9000;
x=1:
16;
w=（x-1）*2*pi/length（x）;
w2=0:
2*pi/256:
255*2*pi/256;
fori=1:
m
t=x/fs（i）;
y=cos（a*t）+cos（b*t）+cos（c*t）;
subplot（6,2,i*2-1）,plot（w,abs（fft（y）））;
y=[yzeros（1,256-16）];
subplot（6,2,i*2）,plot（w2,abs（fft（y）））;
end
共画了6*2个图,每行表示不同的采样率,每行左边表N=16的FFT,右边表示M=256的FFT
结果图如下
对于结果的讨论:
对于序列补0后的FFT与补0前相比,感觉上频谱更"细腻"了点,但实际上,补0后的FFT只是补0前FFT在频域上的插值而已,不会有任何更多的细节出来
另外,当fs=17000,16000时,频谱已不能分辨cos（ct）因为此时的采样频率fs小于2倍cos（c*t）的频率,这与Nyquist采样定理相符合
2、参考以下程序段，对DFT和FFT算法做计算效率比较
Nmax=2048;
fft_time=zeros（1,Nmax）;
forn=1:
Nmax
x=rand（1,n）;
t=clock;fft（x）;fft_time（n）=etime（clock,t）;
end
n=[1:
Nmax];
plot（n,fft_time,'.'）
xlabel（'N'）;ylabel（'TimeinSec.'）
title（'FFTexecutiontimes'）
比如，从N1点到N2点内的效率比较，设N1=8，N2可以考虑取256，或更大些。
大致需要这样一个循环
forn=N1:
N2
产生数据（可以用固定数据，或是随机数据）
计时开始
DFT
计时结束，统计
计时开始
FFT
计时结束，统计
end
显示统计数据
计时函数clock,tic等参阅联机帮助。
此题要自己编写DFT的程序
DFT的函数如下
functionout=dft（x）
m=length（x）;
out=zeros（1,m）;
t=0:
m-1;
t=t*2*pi/m;
forn=1:
m
out（n）=sum（x.*exp（（n-1）*i*t））;
end
end
而主程序段如下
Nmax=256;
ctime=zeros（2,Nmax-63）;
n=1;
forn=64:
Nmax
x=rand（1,n）;
tic;fft（x）;ctime（1,n-63）=toc;
tic;dft（x）;ctime（2,n-63）=toc;
end
n=[64:
Nmax];
plot（n,ctime（1,:
）,'b.'）
holdon,plot（n,ctime（2,:
）,'r-'）
xlabel（'N'）;ylabel（'TimeinSec.'）
legend（'fft','dft'）;
title（'FFTVSDFT'）
结果
3、对抽取/内插前后的信号做傅里叶分析
本次实验的信号均假定起始时间下标为0，也就是对信号x（n）作N:
1抽取，只要
y=x（1:
N:
end）
即可，而1:
N内插则为
y=zeros（1,N*length（x））;
y（[1:
N:
length（y）]）=x;
我们要着重观察的是抽取、内插后频谱的改变。
本实验的数据放在updn.mat文件内，执行语句
loadupdn
要用的数据就会载入数组siga和sigb，如何获取数组大小的信息？
对siga做抽取，sigb做内插实验，试用
N=2,3,4
做抽取，内插，观察它们频谱的变化。
提示：
做谱分析时补一定的0在序列尾部。
实验3FIR数字滤波器的设计
一、实验目的
掌握用MATLAB设计FIR数字滤波器的方法，重点要掌握窗函数法。
二、实验原理
用MATLAB提供的函数，设计FIR数字滤波器
三、实验内容与要求
1.参考示范程序
窗法：
Examples7.8-11（在E:
\DSP-A\RefProgram\CHAP_07目录下）
最优化设计（Parks-McClellan-Remez）：
Examples7.23-25（在E:
\DSP-A\RefProgram\CHAP_07目录下）
在阅读这些例题时，应该注意：
A）如何从滤波器指标要求导出其它设计参数，如确定窗口类型、滤波器的阶数等；
B）例题中验证设计得到的滤波器是否满足设计指标的语句；
C）特别是在最优化设计的例中的迭代设计过程；
D）把滤波器的设计和其特性的显示分开，你自己的显示不一定要照搬书上的。
特别最优化设计的那几个例子程序，显示部分需要函数ampl_res（），需要自己写。
在这基础上，试着直接用我们讲的MATLAB函数fir1（）、fir2（）进行设计
2.用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）观察滤波效果。
这个实验一般不须在时域观察，主要应在频域作观察。
用
freqz（x,1）
就能观察有限长序列的DTFT，对对数显示不习惯的，可以用
[H,w]=freqz（x,1）；
得到数据，再用
freqzplot（H,w,'linear'）
等方式显示。
具体用法参见参考程序。
3．最优化设计方法为选做项目。
四、参考程序：
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,Hb3];
%用不同的方式显示，如果三种都要同时看，你还得加几句
%例如用pause在时间上隔开，或是用3个窗口
figure
（1）
freqzplot（H,w）
figure
（2）
freqzplot（H,w,'linear'）
figure（3）
freqzplot（H,w,'squared'）
end
五、实验所用数值处理函数表
MATLAB中提供了一些数值处理函数，在编程时常会用到，它们也是按照数组运算规则进行运算的。
函数名
功能
函数名
功能
函数名
功能
fix
向零取整
ceil
向正无穷方向取整
mod
除法求余
（与除数同号）
round
四舍五入
floor
向负无穷方向取整
rem
除法求余
（与被除数同号）
我们的网站：
出品：
天涯工作室
联系我们：
2447417264
升级会员 