数字信号处理实验指导书Word下载.docx
《数字信号处理实验指导书Word下载.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书Word下载.docx（36页珍藏版）》请在冰点文库上搜索。
4.我们在这个impseq.m文件中输入如下语句并保存，注意以%开头来写注释。
function[x,n]=impseq（np,ns,nf）
%生成x（n）=delta（n-np）;
ns≤n≤nf，即单位冲激序列δ（n）
%np代表脉冲的位置,ns为序列起始位置，nf为序列终止位置。
%调用方式[x,n]=impseq（np,ns,nf）
ifns>
np|ns>
nf|np>
nf
error（'
输入位置参数不满足ns=<
n<
=nf'
）
elsen=[ns:
nf];
x=[（n-np）==0];
%if判断语句到这里结束
%上面的这条语句是关键语句，作用是如果满足条件（n-np）==0，则x=1。
%在n=ns:
nf的一串值中，只有一个值会满足这个逻辑式，因而只在这个n=np
%处,x为1，其余的n值处，x均为0.这样就构成了延迟np的单位冲激序列
end
5.下面我们再回到图1所示的窗口
1）打算显示一个序列x（n）=1.5*δ（n+1）－δ（n－3）,下面的操作都在命令窗口（CommandWindow）中进行。
注意：
如果一个语句后面有分号，就不显示中间结果。
n1=[-4:
5];
x1=1.5*impseq（-1,-4,5）-impseq（3,-4,5）;
%列出x1序列
subplot（2,2,1）;
stem（n1,x1,'
.'
）;
title（'
x（n）的序列图'
ylabel（'
x1（n）'
axis（[-5,5,-2,3]）;
text（5.5,-2,'
n'
显示的结果如图3所示。
图3序列x（n）=1.5*δ（n+1）－δ（n－3）
2）打算显示另外一个序列
x（n）=n[u（n）－u（n－8）]－10exp（－0.3（n－10））*[u（n－10）—u（n－16）]，0≤n≤20
新建一个文件stepseq.m,用于生成延迟的单位阶跃序列，其输入参数为序列起始位置ns,序列终止位置nf，及阶跃位置np.将刚才编辑的impseq.m另存为stepseq.m，并改动两条语句：
function[x,n]=stepseq（np,ns,nf）
x=[（n-np）>
=0];
存盘。
下面的操作都在命令窗口（CommandWindow）中进行。
n2=[0:
20];
x21=n2.*（stepseq（0,0,20）-stepseq（8,0,20））;
%列出x21序列
x22=10*exp（-0.3*（n2-10））.*（stepseq（10,0,20）-stepseq（16,0,20））;
%列出x22序列
x2=x21-x22;
%x2序列是x21和x22之和
subplot（2,2,2）;
stem（n2,x2,'
x2（n）'
参见附录Matlab简介
实验二离散系统分析
一、实验目的：
1．熟悉和掌握差分方程的求解方法。
2．熟悉和掌握离散付立叶变换和快速付立叶的算法和实现。
3．熟悉和掌握卷积运算。
二、实验内容：
在实验二中，主要有差分方程的求解，DFT及FFT的实现，信号卷积。
差分方程可表示为：
在设计中用到的方程是：
初始条件是：
在实验过程中可以改变系数，在这里只是为了简单，同学们可以在程序中修改一些参数来观察输出的变化。
DFT的实现更简单，在MATLAB语言中有专门的函数，在这里我们选择的输入信号是最简单的余弦函数：
，至于序列的长度的选择，同学们可以在实验中自己选择，建议不要选择太长的序列长度，以免在'
处理过程中影响速度。
建议选择12—64之间。
'
FFT是DFT的快速算法，它的执行速度远远快于DFT，所以可以选择较大的序列长度。
信号卷积为
，在MATLAB中，提供了卷积函数conv，调用十分方便。
在实验中选定的两个信号为：
系统单位冲击响应：
输入信号：
系统的序列长度可以选择的范围是50到100，信号的序列长度选择范围是20到50，采样频率一般选择在1000到3000Hz，至于参数选定，同学们可以自己来决定，实验中用到是
当然在以上设计中信号的选择可以是随意的，并不是固定的。
三、实验报告要求：
1．选择适当数据完成上述实验。
2．给出实验结果并进行分析。
四、参考程序：
参见Untitled.m文件。
实验三用FFT进行信号的频谱分析
1．在理论学习的基础上，通过实验，加深对快速傅立叶变换的理解，熟悉FFT算法；
2．熟悉应用FFT对典型信号进行频谱分析的方法；
3．了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。
在实验中，我们分析三个序列的频谱，分别是高斯序列，衰减正弦序列，还有一般序列。
高斯序列表达式如下：
其中，
可以在实验中自己选定，一般选择规则如下，
选择范围6-10，
为2-8，序列长度选择为16-32，因为太长的话，处理速度会相对比较慢。
衰减正弦序列表示如下：
可以选定的范围为0-1之间，序列长度选择范围为16-32，采样选择为1000-2000之间。
实验中一般序列采用的是：
其中
，
可选的范围为1000-2000,序列长度为16-32。
1）编写程序，实现基二DIT（时域抽取）FFT。
新建一个文件myditfft.m,输入如下语句：
functiony=myditfft（x）
%
%用MATLAB语言编写的基2DITFFT子程序
%
%y=myditfft（x）
%------------------------------------------------------------
%本程序对输入序列x实现DIT-FFT基2算法，
%点数取大于等于x长度的2的幂次
%x为给定时间序列
%y为x的离散傅立叶变换
m=nextpow2（x）;
N=2^m;
%求x的长度对应的2的最低幂次m
iflength（x）<
N
x=[x,zeros（1,N-length（x））];
%若x的长度不是2的幂，补零到2的整数幂
nxd=bin2dec（fliplr（dec2bin（[1:
N]-1,m）））+1;
%求1:
2^m数列的倒序
y=x（nxd）;
%将x倒序排列作为y的初始值
formm=1:
m
%将DFT作m次基2分解,从左到右，对每次分解作DFT运算
Nmr=2^mm;
u=1;
%旋转因子u初始化为WN^0=1
WN=exp（-i*2*pi/Nmr）;
%本次分解的基本DFT因子WN=exp（-i*2*pi/Nmr）
forj=1:
Nmr/2%本次跨越间隔内的各次蝶形运算
fork=j:
Nmr:
N%本次蝶形运算的跨越间隔为Nmr=2^mm
kp=k+Nmr/2;
%确定蝶形运算的对应单元下标
t=y（kp）*u;
%蝶形运算的乘积项
y（kp）=y（k）-t;
%蝶形运算
y（k）=y（k）+t;
end
u=u*WN;
%修改旋转因子,多乘一个基本DFT因子WN
end
[对这个程序的说明]
（1）y=nextpow2（x）:
如果x是一个数，用来求最靠近x并大于x的二的幂次。
例如nextpow2（9.5）=4，因为最靠近9.5，并比它大的2的乘幂是2^4=16;
如果x是一个向量，则nextpow2（x）=nextpow2（length（x））。
即求x长度最近的二的幂次。
（2）Matlab信号处理工具箱中也有倒序函数bitrevorder,为弄清意义，可以自编。
（3）dec2bin（x）把十进制数转换为二进制。
bin2dec（x）把二进制数转换为十进制。
fliplr（x）把x数组排列左右反转。
（4）变量mm的含义相当于FFT分段的序数，最大相当于log2N
2）验证myditfft的正确性，下面的操作都在命令窗口（CommandWindow）中进行。
K=input（'
K='
%设定数据长度的二的幂次K
x=randn（1,2^K）;
%q先生成一个x向量
tic,X=myditfft（x）,toc
%验证myditfft的正确性，并测试所需运行时间。
实验四用窗函数设计FIR数字滤波器
1．熟悉FIR滤波器设计的基本方法；
掌握用窗函数设计FIR数字滤波器的原理及方法；
熟悉线性相位FIR滤波器的幅频特性和相位特性；
2．了解各种不同窗函数对滤波器性能的影响。
这次实验是为了了解各种窗函数的特性，所以实验中用到的窗函数比较多，分别如下：
矩形窗，三角窗，汉宁窗，海明窗，布拉克曼窗，凯泽窗。
关于这几个窗函数的表达形式同学们可以在教材中查到，这里就不一一列举出来了，实验中信号的采样频率可以选择为1000-3000Hz，数字滤波器的截至频率为200Hz，阶数选择范围为41-100。
窗函数选择原则
（1）具有较低的旁瓣幅度，尤其是第一旁瓣的幅度；
（2）旁瓣的幅度下降的速率要快，以利于增加阻带的衰减；
（3）主瓣的宽度要窄，这样可以得到比较窄的过渡带。
实验中还用到了一个特例，具体指标如下：
其中阶数的选择范围为32-100之间
3．简单说明用不同窗函数设计的滤波器之间的性能差异；
实验五设计IIR数字滤波器
1、了解两种工程上最常用的变换方法：
脉冲响应不变法和双线性变换法。
2、掌握双线性变换法设计IIR滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR数字滤波器的计算机程序。
3、观察用双线性变换法设计的滤波器的频域特性，并与脉冲响应不变法相比较，了解双线性变换法的特点。
（i）在各种信号序列中，，以巴特渥斯低通数字滤波器设计为例，可以将双线性变换法设计数字滤波器的步骤归纳如下：
1、确定数字滤波器的性能指标。
这些指标包括：
通带、阻带临界频率
；
通带内的最大衰减
阻带内的最小衰减
采样周期T。
2、确定相应的数字频率
3、计算经过频率预畸的相应参考模拟低通原型的频率，
。
4、计算低通原型阶数N；
计算3dB归一化频率Ωc，从而求得低通原型的传递函数Ha（s）。
5、用变换公式
，代入Ha（s），求得数字滤波器传递函数
.
6、分析滤波器频域特性，检查其指标是否满足要求。
（ii）自编写实现双线性变换的文件bilinear0.m，内容如下：
function[bd,ad]=bilinear0（ba,aa,Fs）
%双线性系数变换子程序
%[bd,ad]=bilinear0（ba,aa,Fs）
%-------------------------------------------------------
%从s域到z域的双线性频带变换
%实现:
%b（z）b（s）|
%----=----|1-z^（-1）Nz（z）
%a（z）a（s）|@s=2Fs--------=------
%1+z^（-1）Dz（z）
%ba,aa按s的正幂降序排列至常数项结束
%Nz,Dz按z的负幂降序排列，从常数项开始，双线性变换：
Nz＝2*Fs*[1,-1];
Dz=[1,1];
%bd,ad按z的负幂降序排列，从常数项开始
%从常数项开始，把ba,aa按s的负幂降序排列。
即将ba,aa中的短者左边补零成同长
lba=length（ba）;
laa=length（aa）;
ld=laa-lba;
ifld>
=0ba=[zeros（1,ld）,ba];
%若aa长度大于ba,給ba前补ld个零
elseaa=[zeros（1,-ld）,aa];
%若ba长度大于aa,給aa前补ld个零
Nz=2*Fs*[1,-1];
Dz=[1,1];
%双线性变换分子分母系数向量
N=max（lba,laa）-1;
%模拟系统阶次N
bd=0;
ad=0;
%bd,ad系数向量初始化
fork=0:
N
pld=[1];
pln=[1];
%双线性变换分子分母系数乘积项初始化
forl=0:
k-1
pld=conv（pld,Dz）;
%求双线性变换分母系数k次幂Dz^k
N-k-1
pln=conv（pln,Nz）;
%求双线性变换分子系数（N-k）次幂Nz^（N-k）
bd=bd+ba（k+1）*conv（pln,pld）;
%分子系数多项式向量求和
ad=ad+aa（k+1）*conv（pln,pld）;
%分母系数多项式向量求和
ad1=ad
（1）;
ad=ad/ad1;
bd=bd/ad1;
%用分母系数多项式首项使分子分母系数归一化
（iii）用巴特沃斯滤波器原型设计一个低通滤波器，满足：
ωp=0.2π，Rp=1dB;
ωs=0.3π，As=15dB;
滤波器的采样频率为1000Hz
自编写实现此功能的文件hc835.m，内容如下：
%双线性变换法由巴特沃思变数字滤波器
%数字滤波器指标:
wp=0.2*pi;
%数字通带频率（弧度）
ws=0.3*pi;
%数字阻带频率（弧度）
Rp=1;
%通带波动（dB）
As=15;
%阻带衰减（dB）
%模拟原型指标的频率逆映射
T=0.001;
Fs=1/T;
%置T=0.001
OmegaP=（2/T）*tan（wp/2）;
%原型通带频率预修正
OmegaS=（2/T）*tan（ws/2）;
%原型阻带频率预修正
ep=sqrt（10^（Rp/10）-1）;
%通带波动参数
Ripple=sqrt（1/（1+ep*ep））;
%通带波动
Attn=1/（10^（As/20））;
%阻带衰减
%模拟巴特沃思原型滤波器计算:
[N,OmegaC]=buttord（OmegaP,OmegaS,Rp,As,'
s'
%原型的阶数和截止频率计算
%%***巴特沃思滤波器阶次=6
[z0,p0,k0]=buttap（N）;
%归一化巴特沃思原型设计函数
p=p0*OmegaC;
z=z0*OmegaC;
%将零极点乘以Omegac，得到非归一化零极点
k=k0*OmegaC^N;
%将k0乘以Omegac^N,得到非归一化k
ba0=real（poly（z0））;
ba0=k0*ba0%由零点计算分子系数向量
aa0=real（poly（p0））%由极点计算分母系数向量
ba=real（poly（z））;
ba=k*ba%由零点计算分子系数向量
aa=real（poly（p））%由极点计算分母系数向量
[bd,ad]=bilinear（ba,aa,Fs）;
%双线性变换:
[bd1,ad1]=bilinear（ba0,aa0,Fs/OmegaC）;
[sos,G]=tf2sos（bd,ad）%变为二阶环节级联结构
%绘图
figure
（1）;
subplot（1,1,1）
[db,mag,pha,grd,w]=myfreqz（bd1,ad1）;
%检验频率响应
plot（w/pi,mag）;
title（'
幅度响应'
xlabel（'
ylabel（'
|H|'
axis（[0,1,0,1.1]）;
set（gca,'
XTickMode'
'
manual'
XTick'
[0,0.2,0.3,1]）;
grid%画刻度线
%set（gca,'
YTickmode'
YTick'
[0,Attn,Ripple,1]）
subplot（2,2,3）;
plot（w/pi,db）;
模值（dB）'
频率：
（单位：
pi）'
分贝'
axis（[0,1,-40,5]）;
%画刻度线
[-50,-15,-1,0]）;
grid
YTickLabelMode'
YTickLabels'
['
50'
;
15'
1'
0'
]）
plot（w/pi,pha/pi）;
相位响应'
单位:
pi'
axis（[0,1,-1,1]）;
[-1,0,1]）;
subplot（2,2,4）;
plot（w/pi,grd）;
群延迟'
频率（单位：
样本'
axis（[0,1,0,10]）
[0:
2:
10]）;
set（gcf,'
color'
w'
）%置图形背景色为白
还需要编写函数myfreqz，新建一个文件myfreqz.m，输入如下内容：
function[db,mag,pha,grd,w]=myfreqz（b,a）;
%freqz子程序的扩展版本，可得出分贝幅特性和群迟延特性
%[db,mag,pha,grd,w]=myfreqz（b,a）;
%------------------------------------
%db=[0到pi弧度]区间内的相对振幅（db）
%mag=[0到pi弧度]区间内的绝对振幅
%pha=[0到pi弧度]区间内的相位响应
%grd=[0到pi弧度]区间内的群延迟
%w=[0到pi弧度]区间内的501个频率样本向量
%b=Ha（z）的分子多项式系数（对FIRb=h）
%a=Ha（z）的分母多项式系数（对FIR:
a=[1]）
[H,w]=freqz（b,a,1000,'
whole'
%用MATLAB中的freqz函数计算数字系统频率响应
H=（H（1:
1:
501））'
w=（w（1:
%取其前一半，并化为列向量
mag=abs（H）;
%求其幅特性
db=20*log10（（mag+eps）/max（mag））;
%化为分贝值
pha=angle（H）;
%求其相特性
grd=grpdelay（b,a,w）;
%求其群迟延特性
程序运行之后，可以观察到频率响应曲线如图所示。
图:
频率响应曲线
实验六随机功率谱估计及MATLAB实现
1．了解功率谱估计的几种算法，比较它们之间的优缺点；
2．熟悉谱估计的算法，以便在实际中应用。
实验中用到的方法有周期性图法，Bartlett平均周期图法，Weleh法及一般方法。
在周期图法中，采样频率
选择范围为1000-3000Hz，FFT算法取样点选择范围为512-1024。
用到的含有噪声的随机序列为：
表示随机序列是Matlab语言自动产生的。
在Bartllett平均周期图法中，采样频率
选择范围为1000-3000Hz，FFT算法取样点选择范围为512-1024。
窗函数阶数选择
升级会员 