数字信号处理实验王翔Word文档格式.docx
- 文档编号:18945717
- 上传时间:2023-01-02
- 格式:DOCX
- 页数:37
- 大小:352.85KB
数字信号处理实验王翔Word文档格式.docx
《数字信号处理实验王翔Word文档格式.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验王翔Word文档格式.docx(37页珍藏版)》请在冰豆网上搜索。
实验地点:
实验结果:
(包括列出实验编写的所有文件及各项实验结果的曲线,并加注必要的说明)
结果分析:
总结:
四、实验考核
(1)实验预习报告;
(2)实验签到;
(3)上机实际操作;
(4)实验设计报告;
五、实验仪器设备要求
(1)有快速的较高性能微机和较大内存与硬盘的设备;
(2)设备数量能适应学生人数;
(3)有Matlab程序设计环境;
六、教材及参考书
1.张志勇等.精通MATLAB6.5.北京航空航天大学出版社.2003年3月
2.邹鲲等.MATLAB6.x信号处理.清华大学出版社.2002年5月
3.陈怀琛等.MATLAB及在电子信息课程中的应用.电子工业出版社.2002
4.程佩清.数字信号处理教程[M].清华大学出版社.2003
实验一离散信号的时域运算与变换
1.熟悉MATLAB编程特点
2.了解离散序列的延迟、相加、相乘及平移、反折、及倒相变换
1.设计一个实现序列移位的函数
将序列x(n)={1,2,3,4,0,7},其中n=0:
5的每一个样本都移动3个周期,移位后的序列y(n)=x(n-3)
2.序列的奇偶分解
将序列x(n)={0,1,2,3,4,3,2,1,0},其中n=-3:
5进行奇偶分解
用函数stem显示其奇偶序列
3.序列的加法运算
设x1(n)={1,0.5,0.3,0.4}其中n=-1:
2;
x2(n)={0.2,0.3,0.4,0.5,0.8,1}其中n=-2:
3
4.实现2中序列的翻转
涉外学院第四教学楼
首先,在用MATLAB表示离散序列并将其可视化时,我们还应注意以下几点:
第一、离散时间序列无法用符号运算来表示,要用矩阵的形式;
第二、由于在MATLAB中矩阵的元素个数有限,因此,MATLAB无法表示无限长的序列;
第三、在绘制离散信号波形的函数stem命令,而不是plot命令。
下面是MATLAB绘制单位抽样序列及其移位序列的函数,n0是单位抽样序列的位移量,n1,n2是序列的起止时刻,,调用该函数就可以绘出单位抽样序列及其移位序列的波形图。
function[x,n]=impuls(n0,n1,n2)
%Generatesx(n)=delta(n-n0);
n=n0处建立一个单位抽样序列
%[x,n]=impuls(n0,n1,n2)
if((n0<
n1)|(n0>
n2)|(n1>
n2))
error('
argumentsmustsatisfyn1<
=n0<
=n2'
)
end
n=[n1:
n2];
x=[zeros(1,(n0-n1)),1,zeros(1,(n2-n0))];
将上述文件存为:
impuls.m,在命令窗口输入
>
n0=1,n1=-20,n2=30;
[x,n]=impuls(n0,n1,n2);
figure
(1),stem(n,x)
上图为右移了一个单位的单位抽样序列
其次,对于离散序列来说,序列相加、相乘是将两序列对应时间序号的值逐项相加或相乘,平移、反折、及倒相变换与连续信号的定义完全相同,但需要注意,与连续信号不同的是,在MATLAB中,离散序列的时域运算和变换不能用符号运算来实现,必须用向量表示的方法,即在MATLAB中离散序列的相加、相乘需表示成两个向量的相加、相乘,因而参加运算的两序列向量必须具有相同的维数。
•在MATLAB中,用一个向量即可表示一个有限长度的序列。
但是,这样的向量并没有包含对应的时间序号的信息。
所以,要完整地表示离散信号要用两个向量
•如:
序列:
f(k)={1,2,-1,3,4,5}
•
k=0
在MATLAB中应表示为:
•f=[1,2,-1,3,4,5]
•k=[-3,-2,-1,0,1,2]或是k=-3:
2
MATLAB常用信号生成函数:
•★ZEROS
•功能:
产生全零阵列
•调用格式:
X=ZEROS(N)%产生N行N列的全零矩阵
•X=ZEROS(M,N)%产生M行N列的全零矩阵
•★
ONES
产生全1阵列
X=ONES(N)%产生N行N列的全1矩阵
•X=ONES(M,N)%产生M行N列的全1矩阵
SINC
辛格函数
Y=SINC(X)%
RECTPULS
产生矩形脉冲信号
Y=RECTPULS(T)%产生高度为1、宽度为1、关于T=0对称的矩形脉冲
•Y=RECTPULS(T,W)%产生高度为1、宽度为W、关于T=0对称的矩形脉冲
RAND
产生伪随机序列
Y=RAND(1,N)%产生[0,1]上均匀分布的随机序列
•Y=RANDN(1,N)%产生均值为0,方差为1的白噪声序列
SAWTOOTH
产生周期锯齿波或三角波
Y=SAWTOOTH(T)%产生幅值为+1,-1,以2为周期的方波
•Y=SAWTOOTH(T,WIDTH)%产生幅值为+1,-1,以WIDTH*2为周期的方波
SQUARE
产生方波
Y=SQUARE(T)%产生幅值为+1,-1,以2为周期的锯齿波
•Y=SQUARE(T,DUTY)%产生幅值为+1,-1,以占空比为DUTY的方波
•例:
t=0:
.0001:
.0625;
•y=SQUARE(2*pi*30*t,80);
plot(t,y)%产生一个占空比为80%的方波
FLIPLR
序列左右翻转
Y=FLIPLR(X)
•%X=123翻转后321
•456654
CUMSUM、SUM
计算序列累加
Y=CUMSUM(X)%向量X元素累加,记录每一次的累加结果,而SUM只记录最后的结果
参考程序:
(1)实现序列移位的函数
function[y,n]=segshift(x,n,n0)
•%功能y(n)=x(n-n0)
•%[y,n]=segshift(x,n,n0)
•n=n+n0;
y=x;
(2)实现两序列相加
function[y,n]=sigadd(x1,n1,x2,n2)
•%实现y(n)=x1(n)+x2(n)
•%-----------------------------
•%[y,n]=sigadd(x1,n1,x2,n2)
•%y=在包含n1和n2的n点上求序列和,
•%x1=在n1上的第一序列
•%x2=在n2上的第二序列(n2可与n1不等)
•n=min(min(n1),min(n2)):
max(max(n1),max(n2));
%y(n)的长度
•y1=zeros(1,length(n));
y2=y1;
%初始化
•y1(find((n>
=min(n1))&
(n<
=max(n1))==1))=x1;
%具有y的长度的x1
•y2(find((n>
=min(n2))&
=max(n2))==1))=x2;
%具有y的长度的x2
•y=y1+y2;
%序列相加.
(3)实现序列奇偶分解
m=-fliplr(n);
•m1=min([m,n]);
m2=max([m,n]);
m=m1:
m2;
•x1=[zeros(1,(length(m)-length(n))),x];
•xe=0.5*(x1+fliplr(x1))
•xo=0.5*(x1-fliplr(x1))
•subplot(1,2,1),stem(m,xe),ylabel('
xe'
),
•subplot(1,2,2),stem(m,xo),ylabel('
xo'
),
奇偶分解结果入下图所示:
4)正弦序列翻转
n=-1:
10;
x=sin(0.4*pi*n);
y=fliplr(x);
n1=-fliplr(n);
subplot(2,1,1),stem(n,x)
subplot(2,1,2),stem(n1,y)
实验二因果离散线性系统的时域分析
实现由差分方程构成的数字滤波器:
1.用两种不同的程序计算数字滤波器的单位抽样响应
,
给定差分方程:
y(n)-y(n-1)+0.9y(n-2)=x(n)
a.编制文件himpz.m,实现数字滤波器的单位抽样响应
b.编制函数文件hfilter.m,实现数字滤波器的单位抽样响应
2.给定
计算此数字滤波器的单位阶跃响应
,并判断系统的稳定性。
用stem(n,y)画出相应的图形。
参考程序gfilter.m
LTI离散系统的单位抽样响应定义为:
当输入为单位抽样序列是系统产生的零状态响应,用h(n)表示。
(1)MATLAB提供了函数impz()求离散系统的单位抽样响应,并绘制时域波形。
★IMPZ
功能:
函数将绘出由矢量a和b表示的数字滤波器在指定取样时刻范围内的抽样响应h(n)的时域波形,并求其数值解。
调用格式:
H=IMPZ(B,A)%计算离散系统冲激响应,取样点数n由函数自动选取
或H=IMPZ(B,A,N)%计算指定时刻范围内(0:
N-1)的离散系统抽样响应
说明:
由矢量a和b构成的数字滤波器的差分方程为
其中:
B=[b0,b1,…,bM,bM-1],A=[a0,a1,…,aN,aN-1],要求A的首项不为0
(2)由于系统的单位冲激响应h(n)也就是系统输入为时系统的零状态响应,除了用上述的impz求解外。
还可以调用filter函数求h(n),此时系统的输入为单位抽样序列。
离散系统差分方程求解的filter函数
调用格式
(1)y=filter(b,a,x)
b,a是式
的系数组成的向量b=[b0,b1,b2,…bM-1]及a=[a0,a1a2,…aN-1],x是输入向量数组,y是输出向量数组和x的长度相同,而且系数a0要保证不为零。
调用格式
(2)y=filter(b,a,x,xic)
xic是初始条件等效的输入序列,MATLAB提供了函数filtic()来确定xic
xic调用格式xic=filtic(b,a,x,y)
b,a是差分方程的系数数组,y,x是初始条件数组,分别由y(n),x(n)的初始条件确定,即y=[y(-1),y(-2),….]及x=[x(-1),x(-2),…]
(3)要判断系统稳定性,一个方法是必须对h(n)求累加,如果结果不为inf(无穷大),则系统稳定;
另一个方法是求系统的极点(即差分方程等号右端为零时方程的解),判断极点是否都在单位园内,即模值是否小于1。
MATLAB提供了一个函数ROOTS(C)求齐次方程的根
★函数ROOTS
求多项式的根
X=ROOTS(C)%计算形如=0的根,
%其中多项式系数矩阵为C,C有N+1个元素,
参考程序
(1)himpz.m
a=[1,-1,0.9];
b=1;
impz(b,a)
title('
单位抽样响应'
(2)hfilter.m
%求h(n)
x=impuls(0,-10,100);
n=[-10:
100];
h=filter(b,a,x);
subplot(2,1,1);
stem(n,h)
axis([-10,100,-1.1,1.1])
title(‘单位冲激响应’);
xlabel(‘n’);
ylabel(‘h(n)’)
(3)gfilter.m
%求单位阶跃响应
x=[zeros(1,10),ones(1,101)];
n=-10:
100;
s=filter(b,a,x);
subplot(2,1,2);
stem(n,s)
axis([-10,100,-.5,2.5])
title(‘单位阶跃响应’);
ylabel(‘s(n)’)
%判断系统稳定
sum(abs(h))
z=roots(a);
magz=abs(z)
程序运行结果:
ans=14.8337
magz=
0.9487
h(n)累加的结果是一个实数,说明系统稳定;
magz说明极点的模小于1,所以系统稳定。
实验三DFT变换的性质及应用
1.实现信号的DFT变换
2.了解DFT应用:
(1)用DFT计算卷积
(2)用DFT对序列进行谱分析
1.用三种不同的DFT程序计算
的傅立叶变换X(k),并比较三种程序的计算机运行时间
步骤:
a.用for循环语句编制函数文件dft1.m,实现循环计算X(k);
b.编写矩阵运算的函数文件dft2.m,实现矩阵计算X(k);
根据定义:
为序列x(n)的DFT
则
c.调用FFT函数直接计算X(K),如程序函数dft3.m
d.分别利用上述三种不同的方式编写的DFT程序计算序列x(n)的DFT变换X(k),并画出幅频和相频特性,并比较3个程序的运行时间。
利用DFT实现两序列的线性卷积运算,并研究DFT的点数与混叠的关系,并用stem(n,y)画出相应的图形。
参考程序dft4.m
3.讨论序列补零及增加数据长度对信号频谱的影响
(1)求出序列x(n)=cos(0.48n)+cos(0.52n)基于有限个样点n=10的频谱
(2)求n=100时,取x(n)的前10个,后90个设为零,得到x(n)的频谱
(3)增加x(n)有效的样点数,取100个样点得到x(n)的频谱参考程序dft5.m
1.dft.m
function[am,pha]=dft1(x)
N=length(x);
w=exp(-j*2*pi/N);
fork=1:
N
sum=0;
forn=1:
sum=sum+x(n)*w^((k-1)*(n-1));
end
am(k)=abs(sum);
pha(k)=angle(sum);
2.dft2.m
function[am,pha]=dft2(x)
n=[0:
N-1];
k=[0:
nk=n’*k;
wnk=w.^(nk);
Xk=x*wnk;
am=abs(Xk);
pha=angle(Xk);
3.dft3.m
function[amfft,phafft]=dft3(x)
Xk=fft(x);
amfft=abs(Xk);
phafft=angle(Xk);
4.dft4.m
%N1+N2-1=23<
32
N=32;
x=[0:
15];
xx=[x,zeros(1,16)];
h=[ones(1,8),zeros(1,24)];
Xk=fft(xx,N);
Hk=fft(h,N);
Yk=Xk.*Hk;
y=ifft(Yk,N);
n=0:
N-1;
stem(n,y);
holdon
%N=N1=16
N1=16;
x1=[0:
h1=[ones(1,8),zeros(1,8)];
Xk1=fft(x1,N1);
Hk1=fft(h1,N1);
Yk1=Xk1.*Hk1;
y1=ifft(Yk1,N1);
n1=0:
N1-1;
stem(n1,y1,’.’,’m’);
5.dft5.m
%x(n)基于10个样点的频谱
figure
(1)
1:
99];
x=cos(0.48*pi*n)+cos(0.52*pi*n);
n1=[0:
9];
y1=x(1:
10);
stem(n1,y1);
title('
signalx(n),0<
=n<
=9'
);
xlabel('
n'
axis([0,10,-2.5,2.5])
Y1=fft(y1);
magY1=abs(Y1(1:
6));
k1=0:
5;
w1=2*pi/10*k1;
stem(w1/pi,magY1);
title(‘10点DFT'
w/pi'
),axis([0,1,0,10])
%在10个样点的基础上添90个零,得到密度高的频谱
figure
(2)
n3=[0:
y3=[x(1:
10)zeros(1,90)];
%添90个零。
得到100个数据
stem(n3,y3);
=9+90zeros'
axis([0,100,-2.5,2.5])
Y3=fft(y3);
magY3=abs(Y3(1:
51));
k3=0:
50;
w3=2*pi/100*k3;
stem(w3/pi,magY3);
100点DFT'
axis([0,1,0,10])
%增加x(n)有效的样点数,取100个样点
figure(3)
stem(n,x);
=99'
X=fft(x);
magX=abs(X(1:
k=0:
w=2*pi/100*k;
stem(w/pi,magX);
100点DFT);
axis([0,1,0,60])
实验四数字低通巴特沃斯滤波器的设计
掌握IIR数字滤波器的设计方法
IIR数字滤波器的设计,主要采用间接法,即:
首先设计出低通模拟滤波器H(S);
进行频率变换,将其转换为高通、带通、带阻滤波器;
再用脉冲响应不变法或双线性变换法从模拟滤波器转换为数字滤波器。
对单极点的N阶H(S)用部分展开式:
冲激响应不变法取H(S)的单阶极点Sk的指数函数
作为H(Z)的极点
双线性变换法是用
代换H(S)中的S得到H(Z),双线性变换法可完全消除频率混叠失真但存在非线性频率失真,而冲激响应不变法存在混叠失真。
在不同的设计阶段MATLAB的信号处理工具箱都给出了相应的滤波器设计函数,这些函数代表了不同类型的逼近函数的滤波器,常用的有巴特沃斯滤波器、切比雪夫滤波器、椭圆滤波器类型。
设计流程图(下图所示为IIR滤波器设计过程及所用到的MATLAB工具箱中的函数)
滤波器系数
A,B
1.设计一个低通巴特沃斯模拟滤波器,绘制滤波器的的幅频响应及零极点图。
指标如下:
通带截止频率:
WP=1000HZ,通带最大衰减:
RP=3dB
阻带截止频率:
Ws=2000HZ,阻带最小衰减:
Rs=40dB
参考程序butter1.m
2.用冲激响应不变法和双线性变换法将一模拟低通滤波器转换为数字低通滤波器
并图释H(S)和H(Z),采样频率Fs=1000Hz
参考程序butter2.m
IIR设计函数
巴特沃斯滤波器为例介绍设计函数的功能和用法,其它设计函数的用法类似。
(1)求最小阶数N的函数buttord
调用格式1:
[N,Wn]=BUTTORD(Wp,Ws,Rp,Rs,'
s'
✧功能:
求出巴特沃斯模拟滤波器的最小阶数N和频率参数Wn
✧说明:
Wp,Ws,为通带截止频率和阻带截止频率如式所示;
Rp,Rs为通带最大衰减Apass(dB)和阻带最小衰减Astop(dB)如式
'
对应模拟滤波器。
(2)模拟低通滤波器设计函数buttap
[Z,P,K]=BUTTAP(N)
按给定的阶数N设计出巴特沃斯模拟低通滤波器
Z,P,K为返回的N阶模拟滤波器的零点、极点和增益系数。
如:
[式3.5.4]
其中k为增益常数,zj和pi为系统函数H(Z)的N个零点和N个极点。
(3)模拟数
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 实验