数字信号处理A实验指导书.docx
《数字信号处理A实验指导书.docx》由会员分享，可在线阅读，更多相关《数字信号处理A实验指导书.docx（30页珍藏版）》请在冰点文库上搜索。
数字信号处理A实验指导书
数字信号处理（A）
实验指导书
柏逢明
吉林城建学院
2012年10月20日
实验一离散时间信号分析
一、实验目的：
加深对常用离散信号的理解，掌握序列运算和序列构造方法。
二、实验原理：
连续信号采样获得离散时间序列，对序列进行运算和构造典型序列过程。
三、实验内容：
信号是信息的载体，它可以传递不同的信息。
离散时间信号是时间变量离散化的信号。
它是对连续时间信号经过抽样而获得的。
经过采样得到的多个离散点的有序排列的信号可称为序列x（n）。
（一）、序列的运算
1.序列相加与相乘：
已知序列
，分别求两序列相加及相乘
%sequence1.m
w0=pi/20;a=1.05;
n1=[-20:
20];n2=[-20:
20];
subplot（2,2,1）,stem（n1,sin（w0.*n1）,'.k'）;axis（[-20,20,-1,3]）;
text（-10,2,'x1（n）'）;
xlabel（'n'）
subplot（2,2,2）,stem（n2,a.^n2,'.k'）;axis（[-20,20,-1,3]）;
text（-10,2,'x2（n）'）;
xlabel（'n'）
n=[min（min（n1）,min（n2））:
max（max（n1）,max（n2））];
x1=zeros（1,length（n））;
x2=zeros（1,length（n））;
x1（[find（（n>=min（n1））&（n<=max（n1）））]）=sin（w0.*n1）;
x2（[find（（n>=min（n2））&（n<=max（n2）））]）=a.^n2;
y1=x1+x2;
y2=x1.*x2;
subplot（2,2,3）;stem（n,y1,'.k'）;axis（[-20,20,-1,3]）;
text（-10,2,'y1（n）'）;
xlabel（'n'）
subplot（2,2,4）;stem（n,y2,'.k'）;axis（[-20,20,-1,3]）;
text（-10,2,'y2（n）'）;
xlabel（'n'）
2.序列移位（延迟）：
设x（n）＝｛1,2,3,4,5,5,4,3,2,1｝且n＝｛0,1,2,3,4,5,6,7,8,9｝，求出序列x（n）及其移位x（n-2）、x（n+2）的图形。
%sequence2.m
n1=[0:
9];
x=[1,2,3,4,5,5,4,3,2,1];
subplot（3,1,1）;stem（n1,x,'.b'）,axis（[-3,12,0,5]）
text（0,3,'x（n）'）
n2=n1-2;
subplot（3,1,2）;stem（n2,x,'.b'）,axis（[-3,12,0,5]）
text（9,3,'x（n+2）'）
n3=n1+2;
subplot（3,1,3）;stem（n3,x,'.b'）,axis（[-3,12,0,5]）
text（0,3,'x（n-2）'）
3.序列反褶（翻转）：
已知序列
，求x（n）的反褶序列。
解：
x（n）的反褶序列为
，设
，a＝1.1，－10≤n≤20。
序列x（n）的反褶序列x（-n）的Matlab程序实现：
%sequence3.m
a=1.1;
n=[-10:
20];x=a.^n;
subplot（2,1,1）;stem（n,x,'.k'）,axis（[-20,20,0,8]）
text（0,5,'x（n）'）
n=fliplr（-n）;x=fliplr（x）;
subplot（2,1,2）;stem（n,x,'.k'）,axis（[-20,20,0,8]）
text（0,5,'x（-n）'）
（二）、典型序列
1.单位冲激序列（单位抽样序列）：
（1）采用Matlab编程实现：
function[x,nx]=delta（n）
x=（n==0）;
nx=n;
（2）采用Matlab的zeros函数实现：
x=zeros（1,N）;
x
（1）=1
2.单位阶跃序列：
（1）采用Matlab编程实现：
functionu=u（n）
u=[n>=0];
nu=n;
（2）采用Matlab的ones函数实现：
x=ones（1,N）;
3.矩形序列：
function[x,nx]=RNn（N,n）;%设N=5，此时R5（n）
x=[（n>=0）&（n<=N-1）];
nx=n;
4.实指数序列：
%sequence4.m
a1=1.1;a2=0.9;a3=-1.1;a4=-0.9;
n=[-5:
15];
x1=（a1.^n）;x2=（a2.^n）;x3=（a3.^n）;x4=（a4.^n）;
subplot（2,2,1）;stem（n,x1,'.k'）;axis（[-5,15,-0.5,5]）;
xlabel（'n'）;
text（3.5,3.5,'a>1'）;
subplot（2,2,2）;stem（n,x2,'.k'）;axis（[-5,15,-0.2,1.2]）;
xlabel（'n'）;
text（10,0.7,'0subplot（2,2,3）;stem（n,x3,'.k'）;axis（[-5,15,-6,4]）;
xlabel（'n'）;
text（4,-4,'a<-1'）;
subplot（2,2,4）;stem（n,x4,'.k'）;axis（[-5,15,-1,1.2]）;
xlabel（'n'）;
text（10,0.7,'-15.正弦序列：
，设
%sequence5.m
w=pi/15;%正弦序列的数字域频率反映了序列变化的快慢，如w=0.01*pi和w=0.1*pi
n=[-15:
15];
x=sin（w.*n）;
stem（n,x,'.k'）;axis（[-15,15,-2,2]）;
text（-14,1.7,'sin（nw）'）;
text（12,-0.3,'n'）;
6.复指数序列：
设
，复指数序列
幅度特性与相位特性的Matlab实现程序：
%sequence6.m
w0=pi/15;
n=[-15:
15];
x=exp（j.*w0.*n）;
subplot（2,1,1）;stem（n,abs（x）,'.k'）;axis（[-15,15,-0.2,1.5]）;
text（11,-0.1,'n'）;
text（1,1.2,'x（n）'）
subplot（2,1,2）;stem（n,angle（x）,'.k'）;
text（11,-0.5,'n'）;
text（1,3,'arg[x（n）]'）;
四、实验要求：
①要求学生运用MATLAB编程完成一些数字信号处理的基本功能，加深对教学内容的理解；
②讨论复指数序列的性质。
实验二DFS、DFT及FFT
1.M脚本文件编写；2.Fourier级数与FFT
一、实验目的：
熟悉MATLAB软件中M文件的编写方法，实验数字信号处理从Fouriere级数展开到FFT的过程。
二、实验原理：
时域函数的Fourier级数展开与FFT，符号方法与几分方法和子程序调用。
三、实验内容：
（一）、关于Matlab程序设计之M文件编写：
Matlab语言代码的文件称为M文件，其后缀都为m。
M文件分为函数程序和脚本程序两种：
函数程序可以接受输入参数，并可以输出；脚本程序是Matlab表达式的集合，不可以接受参数。
编辑M文件可以使用各种文本编辑器，Matlab中具有内置的M文件编辑器。
1.M脚本文件
脚本程序是一种最简单的M文件，它没有输入输出参数，只不过是Matlab表达式的集合（如前命令行中的所有命令）。
脚本程序可以是一系列命令行中执行命令的集合，也可以是操作工作空间中的变量和程序中新建的变量。
编写M文件的步骤：
1点击File→NEW选M-file或在窗口工具条上点击NewFile图标，打开Matlab文件编辑器（MatlabEditor/Debugger），进行M文件程序编写（当然，也可以用其他文本编辑器）；
2单击编辑器工具条保存图标，进行文件保存为*.m格式文件,如ex91.m；
3运行文件：
a.必须使得*.m文件所在目录为当前目录，或让该目录处于搜索路径上；b.然后运行该指令：
ex91（即扩展名前的名称）。
2.M函数程序
函数程序既可以接受输入参数，也可以输出参数，同时也可进行函数操作工作空间的变量。
同其他高级语言类似，M函数文件也有局部变量和全局变量。
编写M函数程序文件与编写M脚本文件的步骤相同。
注意：
在一个M函数程序文件中要有函数项语句。
只能有一个主函数（function），一般放在文件最前面，其它函数均为子函数（function）。
运行时运行主函数所制定的文件。
3.应用举例
问题：
已知一分段函数
，试画出此函数所表示的三维曲面。
（1）用M脚本文件表示：
%[ex91.m]
a=2;b=2;
clf;
x=-a:
0.2:
a;y=-b:
0.2:
b;
fori=1:
length（y）
forj=1:
length（x）
ifx（j）+y（i）>1
z（i,j）=0.5457*exp（-0.75*y（i）^2-3.75*x（j）^2-1.5*x（j））;
elseifx（j）+y（i）<=-1
z（i,j）=0.5457*exp（-0.75*y（i）^2-3.75*x（j）^2+1.5*x（j））;
elsez（i,j）=0.7575*exp（-y（i）^2-6.*x（j）^2）;
end
end
end
axis（[-a,a,-b,b,min（min（z））,max（max（z））]）;
colormap（flipud（winter））;surf（x,y,z）;
运行命令：
ex91
运行结果：
三维图形
（2）用M程序文件表示：
%[ex92.m]
functionex92（a,b）％这是主函数
%Thisismysecondexample.
%adefinethelimitofvariablex.
%bDefinethelimitofvariabley.
a=2;b=2;
clf;
x=-a:
0.2:
a;y=-b:
0.2:
b;
fori=1:
length（y）
forj=1:
length（x）
ifx（j）+y（i）>1
z（i,j）=0.5457*exp（-0.75*y（i）^2-3.75*x（j）^2-1.5*x（j））;
elseifx（j）+y（i）<=-1
z（i,j）=0.5457*exp（-0.75*y（i）^2-3.75*x（j）^2+1.5*x（j））;
elsez（i,j）=0.7575*exp（-y（i）^2-6.*x（j）^2）;
end
end
end
axis（[-a,a,-b,b,min（min（z））,max（max（z））]）;
colormap（flipud（winter））;surf（x,y,z）;
运行命令：
ex92（2,2）
运行结果：
三维图形（与前相同）
（二）、数据信号处理：
从Fouriere级数展开到FFT：
已知时域函数
，求Fouriere级数展开系数。
编程方法：
（1）符号方法：
%
%[fzzysym.m]
function[A_sym,B_sym]=fzzysym（T,Nf,Nn）
%采用符号计算求[0,T]内时间函数的三角级数展开系数
%函数的输入输出都是数值量
%Nf谐波的阶数
%Nn输出数据的准确位数
%A_sym第1元素是直流项，其后元素依次是1,2,3,...次谐波cos项展开系数
%B_sym第2,3,4,...元素依次是1,2,3,...次谐波sin项展开系数
symstttn
ifnargin<2;Nf=6;end
ifnargin<3;Nn=32;end
yy=time_fun_s（ttt）;
A0=int（yy,ttt,0,T）/T;
As=int（yy*cos（2*pi*n*ttt/T）,ttt,0,T）;
Bs=int（yy*sin（2*pi*n*ttt/T）,ttt,0,T）;
A_sym
（1）=double（vpa（A0,Nn））;
fork=1:
Nf
A_sym（k+1）=double（vpa（subs（As,n,k）,Nn））;
B_sym（k+1）=double（vpa（subs（Bs,n,k）,Nn））;
end
functionyy=time_fun_s（ttt）
%该函数是fzzysym.m的子函数。
它由符号变量和表达式写成。
y1=sym（'Heaviside（ttt-0.5）'）*（ttt-0.5）;
yy=y1-sym（'Heaviside（ttt-1.5）'）*（（ttt-1.5）+1）;
%运行：
在命令窗口（CommandWindow）中运行函数命令：
[A_sym,B_sym]=fzzysym
（2）
即可得到Fourier级数展开系数结果。
（2）积分方法：
%[fzzyquad.m]
function[t,y,S,A,B]=fzzyquad（a,T,Nf,K）
%用数值积分法求Fourier级数展开系数A、B和近似波形S。
如调用本命令时，没有任何输出总量，那么将自动随积分进程动态地画出原波形和各阶近似波形。
%
%a是被展开函数的时间区间的左端。
%T是被展开函数的周期。
%Nf是所需展开的最高谐波阶次
%K绘制波形所采用的数据长度
%t,y是被展开的时间函数数据
%S是一个矩阵。
它每行的长度与t,y相同。
%它的第（n+1）行是所有n次及更低次谐波叠加生成的近似波形。
%A,B是展开系数向量。
A（n+1）,B（n+1）分别存放cos，sin项n次谐波的展开系数。
%A
（1）是直流项
ifnargin<4;K=100;end
if（nargin<3|isempty（Nf））;Nf=15;end
k=1:
K;
t=a+k*T/length（k）;
y=time_fun（t,T）;%调用函数文件，产生被展开函数。
S=zeros（Nf+1,length（k））;
a0=mean（y）;
S（1,:
）=a0;
A=zeros（1,Nf+1）;B=zeros（1,Nf+1）;
A
（1）=a0;
forn=1:
Nf
A（n+1）=quad8（'cos_y',a,a+T,[],[],n,T）/T*2;
B（n+1）=quad8（'sin_y',a,a+T,[],[],n,T）/T*2;
S（n+1,:
）=S（n,:
）+A（n+1）*cos（2*n*pi*t/T）+B（n+1）*sin（2*n*pi*t/T）;
ifnargout==0
plot（t,y,'b'）;holdon
plot（t,S（n+1,:
）,':
r'）;
aa=axis;dx=aa
（2）-aa
（1）;dy=aa（4）-aa（3）;
text（aa
（1）+0.1*dx,aa（3）+.75*dy,'n='）
text（aa
（1）+0.2*dx,aa（3）+.75*dy,num2str（n））;
pause
（1）;holdoff
end
end
%
%[cos_y.m]余弦项子程序
functionwcos=cos_y（t,n,T）
%生成求cos项展开系数的被积函数。
y=time_fun（t,T）;wcos=cos（2*n*pi*t/T）.*y;
％
%[sin_y.m]正弦项子程序
functionwsin=sin_y（t,n,T）
%生成求sin项展开系数的被积函数。
y=time_fun（t,T）;wsin=sin（2*n*pi*t/T）.*y;
%
%[time_fun.m]时间函数子程序
functiony=time_fun（t,T）
%这是用户自己编写的待展开的时间函数：
y（t+k*T）=y（t）
%编写程序时，一定要注意使表达式适于数组运算
%t是时间数组
%T是周期
y=zeros（size（t））;ii=find（t>=0.5&t<=1.5）;
y（ii）=ones（size（ii））.*（t（ii）-0.5）;y（y==1.0）=0.5;
%注意：
首先在命令窗口中，运行第一行<1>指令：
%[t,y,S,A_quad,B_quad]=fzzyquad（0,2,15）;%计算0到15次谐波<1>
%然后再在命令窗口中运行下面的指令：
%AA=abs（A_sym）;AA（AA<1e-10）=NaN;
%在那些实际应为0的很小数处，考虑相对误差没有意义，故予以排除。
%BB=abs（B_sym）;BB（BB<1e-10）=NaN;
%asq=abs（A_quad（1:
7）-A_sym）./AA
%偶函数系数与符号计算的相对误差
%bsq=abs（B_quad（1:
7）-B_sym）./BB
%奇函数系数与符号计算值的相对误差
%图形显示截断余项后三角展开的近似波形：
%SS=[S;y]';ribbon（t,SS）;%把原始波形列在图形最前方。
得到原始波形。
%view（[96,36]）,colormap（jet）,shadingflat,light,lightinggouraud
%截断后得到的近似波形。
%
（3）FFT方法：
%
function[A,B,C,fn,t,w]=fzzyfft（T,M,Nf）
%利用FFT，计算[0,T]区间上定义的时间波形的Fourier级数展开数A、B和频谱C、Nf。
%T时间波形周期
%M用作2的幂次
%Nf输出谐波的阶次，决定A、B的长度为[Nf+1]。
Nf不要超过2^（M-1）。
%A,B分别是Fourier级数中cos，sin展开项的系数。
A
（1）是直流量。
%C是定义在[-fs/2,fs/2]上的频谱
%t,w是原时间波形数据对
if（nargin<2|isempty（M））;M=8;end%
ifnargin<3;Nf=6;end%
N=2^M;%使总采样点是2的整数倍
f=1/T;%被变换函数的频率
w0=2*pi*f;
dt=T/N;%时间分辨率
n=0:
1:
（N-1）;%采样点序列
t=n*dt;%采样时间序列
w=time_fun（t,T）;%被变换时间函数的采样序列<17>
W=fft（w）;%给出n=0,1,2,....,N-1上的DFT数据值
cn=W/N;%根据公式计算n=0,1,...,N-1上的FS系数
z_cn=find（abs（cn）<1.0e-10）;%寻找有限字长运算而产生（原应为0）的“小”复数
cn（z_cn）=zeros（length（z_cn）,1）;%强制那些“小”复数为0
cn_SH=fftshift（cn）;%根据式（5.13.2.2-3）计算
%n=-N/2,...,-1,0,1,...,（N/2）-1上的FS系数
C=[cn_SHcn_SH
（1）];%形成关于0对称的（N+1）个FS系数
A
（1）=C（N/2+1）;%
A（2:
N/2+1）=2*real（C（（N/2+2）:
end））;
B（2:
N/2+1）=-2*imag（C（（N/2+2）:
end））;
ifNf>N/2;error（['第三输入总量Nf应小于'int2str（N/2-1）]）;end
A（Nf+2:
end）=[];
B（Nf+2:
end）=[];
n1=-N/2:
1:
N/2;%产生总点数为（N-1）关于0对称的序列
fn=n1*f;%关于0对称的频率分度序列
%注意：
在命令窗口运行如下命令：
（一定要首先运行以前的几个相关命令才能接续下来！
）
%[A_fft,B_fft]=fzzyfft
（2）;主程序命令
%AA=abs（A_sym）;AA（AA<1e-10）=NaN;BB=abs（B_sym）;BB（BB<1e-10）=NaN;
%ast=abs（A_sym-A_fft）./AA%偶函数系数与符号计算值的相对误差
%bst=abs（B_sym-B_fft）./BB%奇函数系数与符号计算值的相对误差
四、实验要求：
利用MATLAB编程完成计算，绘出相应图形。
并与理论计算相比较，说明实验结果的原因。
实验三离散系统响应与矩形脉冲时域频谱
一、实验目的：
掌握系统响应和系统特性的实现原理和方法，了解采样频率参数变化对混迭现象的影响。
独立分析矩形脉冲时域频谱。
二、实验原理：
信号系统响应原理，系统特性中零极点和系统频率响应的幅度特性和相位特性定义。
利用DFT计算一般连续函数的傅立叶变换CFT。
三、实验内容：
1.系统响应：
已知离散系统的差分方程为
，求以下集中情况时系统的响应y（n）。
①
；②
；
③
x（n）
设a＝1.01,b=1.1，求解系统响应y（n）的实现程序如下：
%response1.m
%Systemresponse1
clear
a=1.01;
b=1.1;
B=[1,0];
A=[1,-b];
Y=[5];X=[];
n=[0:
30];
x1n=a.^n;
y1n=filter（B,A,x1n）;
subplot（3,1,1）
stem（n,y1n,'.k'）;axis（[0,30,0,300]）
text（10,200,'y1（n）'）
xlabel（'n'）
xic=filtic（B,A,Y,X）;
x2n=zeros（1,length（n））;
y2n=filter（B,A,x2n,xic）;
subplot（3,1,2）
stem（n,y2n,'.k'）;axis（[0,30,0,300]）
text（10,200,'y2（n）'）
xlabel（'n'）
x3n=x1n;
y3n=filter（B,A,x3n,xic）;
subplot（3,1,3）
stem（n,y3n,'.k'）;axis（[0,30,0,300]）
xlabel（'n'）
text（10,200,'y3（n）'）
2.系统特性
已知一阶系统的差分方程为
，求当a=0.9时该系统的频率特性。
设a=0.9,求解系统零、极点的Matlab程序：
clear
b=[1,0];
a=[1,-0.9];
zplane（b,a）;
设a=0.9,求解系统频率响应的幅度特性和相位特性的Matlab程序：
clear
b=[1,0];
a=[1,-0.9];
w=pi*freqspace（100）;
freqz（b,a,w）;
3．采样频率与混迭现象
升级会员 