数字信号处理Matlab课后实验吴镇扬资料.docx

数字信号处理Matlab课后实验吴镇扬资料.docx

《数字信号处理Matlab课后实验吴镇扬资料.docx》由会员分享，可在线阅读，更多相关《数字信号处理Matlab课后实验吴镇扬资料.docx（48页珍藏版）》请在冰豆网上搜索。

数字信号处理Matlab课后实验吴镇扬资料

JIANGSUUNIVERSITY

数字信号处理实验报告

实验一熟悉MATLAB环境

实验二快速变换及其应用

实验三IIR数字滤波器的设计

实验四FIR数字滤波器的设计

实验八信号的谱分析及分段卷

实验一熟悉MATLAB环境

一、实验目的

（1）熟悉MATLAB的主要操作命令。

（2）学会简单的矩阵输入和数据读写。

（3）掌握简单的绘图命令。

（4）用MATLAB编程并学会创建函数。

（5）观察离散系统的频率响应。

二、实验内容

认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。

在熟悉了MATLAB基本命令的基础上，完成以下实验。

上机实验内容：

（1）数组的加、减、乘、除和乘方运算。

输入A=[1234]，B=[3456]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。

实验程序：

A=[1234];

B=[3456];

n=1:

4;

C=A+B;D=A-B;E=A.*B;F=A./B;G=A.^B;

subplot（4,2,1）;stem（n,A,'fill'）;xlabel（'时间序列n'）;ylabel（'A'）;

subplot（4,2,2）;stem（n,B,'fill'）;xlabel（'时间序列n'）;ylabel（'B'）;

subplot（4,2,3）;stem（n,C,'fill'）;xlabel（'时间序列n'）;ylabel（'A+B'）;

subplot（4,2,4）;stem（n,D,'fill'）;xlabel（'时间序列n'）;ylabel（'A-B'）;

subplot（4,2,5）;stem（n,E,'fill'）;xlabel（'时间序列n'）;ylabel（'A.*B'）;

subplot（4,2,6）;stem（n,F,'fill'）;xlabel（'时间序列n'）;ylabel（'A./B'）;

subplot（4,2,7）;stem（n,G,'fill'）;xlabel（'时间序列n'）;ylabel（'A.^B'）;

运行结果：

（2）用MATLAB实现以下序列。

a）x（n）=0.8n0≤n≤15

实验程序：

n=0:

15;x=0.8.^n;

stem（n,x,'fill'）;xlabel（'时间序列n'）;ylabel（'x（n）=0.8^n'）;

b）x（n）=e（0.2+3j）n0≤n≤15

实验程序：

n=0:

15;x=exp（（0.2+3*j）*n）;

stem（n,x,'fill'）;xlabel（'时间序列n'）;ylabel（'x（n）=exp（（0.2+3*j）*n）'）;

运行结果：

a）的时间序列b）的时间序列

c）x（n）=3cos（0.125πn+0.2π）+2sin（0.25πn+0.1π）0≤n≤15

实验程序：

n=0:

1:

15;

x=3*cos（0.125*pi*n+0.2*pi）+2*sin（0.25*pi*n+0.1*pi）;

stem（n,x,'fill'）;xlabel（'时间序列n'）;

ylabel（'x（n）=3*cos（0.125*pi*n+0.2*pi）+2*sin（0.25*pi*n+0.1*pi）'）;

运行结果：

d）将c）中的x（n）扩展为以16为周期的函数x16（n）=x（n+16）,绘出四个周期

实验程序：

n=0:

1:

63;

x=3*cos（0.125*pi*rem（n,16）+0.2*pi）+2*sin（0.25*pi*rem（n,16）+0.1*pi）;

stem（n,x,'fill'）;xlabel（'时间序列n'）;ylabel（'x16（n）'）;

e）将c）中的x（n）扩展为以10为周期的函数x10（n）=x（n+10）,绘出四个周期

实验程序：

n=0:

1:

39;

x=3*cos（0.125*pi*rem（n,10）+0.2*pi）+2*sin（0.25*pi*rem（n,10）+0.1*pi）;

stem（n,x,'fill'）;xlabel（'时间序列n'）;ylabel（'x10（n）'）;

运行结果：

d）的时间序列e）的时间序列

（3）x（n）=[1,-1,3,5]，产生并绘出下列序列的样本。

a）x1（n）=2x（n+2）-x（n-1）-2x（n）

实验程序：

n=0:

3;

x=[1-135];

x1=circshift（x,[0-2]）;x2=circshift（x,[01]）;x3=2*x1-x2-2*x;

stem（x3,'fill'）;xlabel（'时间序列n'）;ylabel（'x1（n）=2x（n+2）-x（n-1）-2x（n）'）;

b）

实验程序：

n=0:

3;

x=[1-135];

x1=circshift（x,[01]）;x2=circshift（x,[02]）;x3=circshift（x,[03]）;

x4=circshift（x,[04]）;x5=circshift（x,[05]）;

xn=1*x1+2*x2+3*x3+4*x4+5*x5;

stem（xn,'fill'）;xlabel（'时间序列n'）;

ylabel（'x2（n）=x（n-1）+2x（n-2）+3x（n-3）+4x（n-4）+5x（n-5）'）;

运行结果：

a）的时间序列b）的时间序列

（4）绘出时间函数的图形，对x轴、y轴图形上方均须加上适当的标注。

a）x（t）=sin（2πt）0≤t≤10sb）x（t）=cos（100πt）sin（πt）0≤t≤4s

实验程序：

clc;

t1=0:

0.001:

10;t2=0:

0.01:

4;

xa=sin（2*pi*t1）;xb=cos（100*pi*t2）.*sin（pi*t2）;

subplot（2,1,1）;

plot（t1,xa）;xlabel（'t'）;ylabel（'x（t）'）;title（'x（t）=sin（2*pi*t）'）;

subplot（2,1,2）;

plot（t2,xb）;xlabel（'t'）;ylabel（'x（t）'）;title（'x（t）=cos（100*pi*t2）.*sin（pi*t2）'）;

运行结果：

（5）编写函数stepshift（n0,n1,n2）实现u（n-n0）,n1

实验程序：

clc;

n1=input（'请输入起点：

'）;

n2=input（'请输入终点'）;

n0=input（'请输入阶跃位置'）;

n=n1:

n2;

x=[n-n0>=0];

stem（n,x,'fill'）;xlable（'时间序列n'）;ylable（'u（n-n0）'）;

请输入起点：

2

请输入终点：

8

请输入阶跃位置：

6

运行结果：

（5）运行结果（6）运行结果

（6）给一定因果系统

求出并绘制H（z）的幅频响应与相频响应。

实验程序：

a=[1-0.670.9];

b=[1sqrt

（2）1];

[hw]=freqz（b,a）;

fp=20*log（abs（h））;

subplot（2,1,1）;

plot（w,fp）;xlabel（'时间序列t'）;ylabel（'幅频特性'）;

xp=angle（h）;

subplot（2,1,2）;

plot（w,xp）;xlabel（'时间序列t'）;ylabel（'相频特性'）;

运行结果：

（右上图）

实验二快速变换及其应用

 一、实验目的

（1）在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉FFT子程序。

（2）熟悉应用FFT对典型信号进行频谱分析的方法。

（3）了解应用FFT进行信号频谱分析过程中可能出现的问题以便在实际中正确应用FFT。

（4）熟悉应用FFT实现两个序列的线性卷积的方法。

（5）对DCT变换用作数据压缩有初步的认识。

二、实验内容

（1）观察高斯序列的时域和幅频特性。

实验程序：

clc;

n=0:

1:

15;

%p=8,q=2

p=8;q=2;

x1=exp（-（n-p）.^2./q）;

fp1=fft（x1）;fp1=abs（fp1）;

subplot（5,2,1）;plot（x1）;xlabel（'n'）;ylabel（'时域特性'）;title（'p=8,q=2'）;

subplot（5,2,2）;plot（fp1）;xlabel（'n'）;ylabel（'幅频特性'）;title（'p=8,q=2'）;

%p=8,q=4

p=8;q=4;

x2=exp（-（n-p）.^2./q）;

fp2=fft（x2）;fp2=abs（fp2）;

subplot（5,2,3）;plot（x2）;xlabel（'n'）;ylabel（'时域特性'）;title（'p=8,q=4'）;

subplot（5,2,4）;plot（fp2）;xlabel（'n'）;ylabel（'幅频特性'）;title（'p=8,q=4'）;

%p=8,q=8

p=8;q=8;

x3=exp（-（n-p）.^2./q）;

fp3=fft（x3）;fp3=abs（fp3）;

subplot（5,2,5）;plot（x3）;xlabel（'n'）;ylabel（'时域特性'）;title（'p=8,q=8'）;

subplot（5,2,6）;plot（fp3）;xlabel（'n'）;ylabel（'幅频特性'）;title（'p=8,q=8'）;

%p=13,q=8

p=13,q=8;

x4=exp（-（n-p）.^2./q）;

fp4=fft（x4）;fp4=abs（fp4）;

subplot（5,2,7）;plot（x4）;xlabel（'n'）;ylabel（'时域特性'）;title（'p=13,q=8'）;

subplot（5,2,8）;plot（fp4）;xlabel（'n'）;ylabel（'幅频特性'）;title（'p=13,q=8'）;

%p=14,q=8

p=14,q=8;

x5=exp（-（n-p）.^2./q）;

fp5=fft（x5）;fp5=abs（fp5）;

subplot（5,2,9）;plot（x5）;xlabel（'n'）;ylabel（'时域特性'）;title（'p=14,q=8'）;

subplot（5,2,10）;plot（fp5）;xlabel（'n'）;ylabel（'幅频特性'）;title（'p=14,q=8'）;

运行结果：

（2）观察衰减正弦序列xb（n）的时域和幅频特性。

=0.1,f=0.0625,检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625,观察这两种情况下频谱的形状和谱峰出现为止，又无混叠和泄露现象？

说明产生现象的原因。

实验程序：

clc;

n=0:

1:

15;

a=0.1;f1=0.0625;f2=0.4375;f3=0.5625;

x1=exp（-a*n）.*sin（2*pi*f1.*n）;

x2=exp（-a*n）.*sin（2*pi*f2.*n）;

x3=exp（-a*n）.*sin（2*pi*f3.*n）;

fp1=fft（x1）;fp1=abs（fp1）;

fp2=fft（x2）;fp2=abs（fp2）;

fp3=fft（x3）;fp3=abs（fp3）;

subplot（3,2,1）;plot（x1）;xlabel（'n'）;ylabel（'时域特性'）;title（'a=0.1,f=0.0625'）;

subplot（3,2,2）;plot（fp1）;xlabel（'n'）;ylabel（'幅频特性'）;title（'a=0.1,f=0.0625'）;

subplot（3,2,3）;plot（x2）;xlabel（'n'）;ylabel（'时域特性'）;title（'a=0.1,f=0.4375'）;

subplot（3,2,4）;plot（fp2）;xlabel（'n'）;ylabel（'幅频特性'）;title（'a=0.1,f=0.4375'）;

subplot（3,2,5）;plot（x3）;xlabel（'n'）;ylabel（'时域特性'）;title（'a=0.1,f=0.5625'）;

subplot（3,2,6）;plot（fp3）;xlabel（'n'）;ylabel（'幅频特性'）;title（'a=0.1,f=0.5625'）;

运行结果：

（3）观察三角波和反三角波序列的时域和幅频特性。

实验程序：

N=8时：

clc;

n=0:

1:

7;

%x1为三角波序列时域特性，x2为反三角波时域特性

x1=zeros（size（n））;x2=zeros（size（n））;

xa=（n>=0）&（n<=3）;x1（xa）=n（xa）;x2（xa）=4-n（xa）;

xa=（n>=4）&（n<=7）;x1（xa）=8-n（xa）;x2（xa）=n（xa）-4;

fp1=fft（x1）;fp1=abs（fp1）;

fp2=fft（x2）;fp2=abs（fp2）;

%fp1为三角波序列幅频特性，fp2为反三角波幅频特性

subplot（2,2,1）;stem（n,x1,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'三角波序列'）;

subplot（2,2,2）;stem（n,fp1,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'三角波序列'）;

subplot（2,2,3）;stem（n,x2,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'反三角波序列'）;

subplot（2,2,4）;stem（n,fp2,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'反三角波序列'）;

运行结果：

实验程序：

N=16时：

clc;

n=0:

1:

16;

%x1为三角波序列时域特性，x2为反三角波时域特性

x1=zeros（size（n））;x2=zeros（size（n））;

xa=（n>=0）&（n<=3）;x1（xa）=n（xa）;x2（xa）=4-n（xa）;

xa=（n>=4）&（n<=7）;x1（xa）=8-n（xa）;x2（xa）=n（xa）-4;

xa=（n>=8）&（n<=15）;x1（xa）=0;x2（xa）=0;

fp1=fft（x1）;fp1=abs（fp1）;

fp2=fft（x2）;fp2=abs（fp2）;

%fp1为三角波序列幅频特性，fp2为反三角波幅频特性

subplot（2,2,1）;stem（n,x1,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'三角波序列'）;

subplot（2,2,2）;stem（n,fp1,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'三角波序列'）;

subplot（2,2,3）;stem（n,x2,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'反三角波序列'）;

subplot（2,2,4）;stem（n,fp2,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'反三角波序列'）;

运行结果：

（4）一个连续信号含两个频率分量，经采样得x（n）=sin2π*0.125n+cos2π*（0.125+Δf）nn=0,1……,N-1已知N=16，Δf分别为1/16和1/64，观察其频谱；当N=128时，Δf不变，其结果有何不同。

实验程序：

N=16时：

clc;n1=0:

1:

15;

%N=16时，vf1=1/16,vf2=1/64两种情况下

vf1=1/16;vf2=1/64;

x1=sin（2*pi*0.125.*n1）+cos（2*pi*（0.125+vf1）.*n1）;%时域特性

x2=sin（2*pi*0.125.*n1）+cos（2*pi*（0.125+vf2）.*n1）;

fp1=fft（x1）;fp1=abs（fp1）;

fp2=fft（x2）;fp2=abs（fp2）;

subplot（2,2,1）;

stem（n1,x1,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'N=16vf=1/16'）;

subplot（2,2,2）;

stem（n1,fp1,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'N=16vf=1/16'）;

subplot（2,2,3）;

stem（n1,x2,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'N=16vf=1/64'）;

subplot（2,2,4）;

stem（n1,fp2,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'N=16vf=1/64'）;

运行结果：

N=128时：

clc;

n2=0:

1:

127;

%N=32时，vf1=1/16,vf2=1/64两种情况下

vf1=1/16;vf2=1/64;

x3=sin（2*pi*0.125.*n2）+cos（2*pi*（0.125+vf1）.*n2）;%时域特性

x4=sin（2*pi*0.125.*n2）+cos（2*pi*（0.125+vf2）.*n2）;

fp3=fft（x3）;fp3=abs（fp3）;%幅频特性

fp4=fft（x4）;fp4=abs（fp4）;

subplot（2,2,1）;

stem（n2,x3,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'N=128vf=1/16'）;

subplot（2,2,2）;

stem（n2,fp3,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'N=128vf=1/16'）;

subplot（2,2,3）;

stem（n2,x4,'.'）;xlabel（'n'）;ylabel（'幅频特性'）;title（'N=128vf=1/64'）;

subplot（2,2,4）;

stem（n2,fp4,'.'）;xlabel（'n'）;ylabel（'时域特性'）;title（'N=128vf=1/64'）;

运行结果：

（5）用FFT分别计算xa（n）（p=8,q=2）和xb（n）（a=0.1,f=0.0625）的16点循环卷积和线性卷积。

实验程序：

clc;

n1=0:

1:

15;

n2=0:

1:

30;

p=8;q=2;a=0.1;f=0.0625;

xa=exp（-（（n1-p）.^2）/q）;

xb=exp（-a.*n1）.*sin（2*pi*f.*n1）;

fa=fft（xa）;fb=fft（xb）;

circle=fa.*fb;%圆周卷积

line=conv（xa,xb）;%线性卷积

subplot（2,2,1）;

stem（n1,xa,'.'）;ylabel（'时域特性'）;title（'高斯序列xa'）;

subplot（2,2,2）;

stem（n1,xb,'.'）;ylabel（'时域特性'）;title（'衰减正弦序列xb'）;

subplot（2,2,3）;

stem（n1,circle,'.'）;ylabel（'幅频特性'）;title（'xa与xb的16点循环卷积'）;

subplot（2,2,4）;

stem（n2,line,'.'）;ylabel（'幅频特性'）;title（'xa与xb的线性卷积'）;

运行结果：

（6）产生一512点的随机序列xe（n），并用xc（n）和xe（n）做线性卷积，观察卷积前后xe（n）频谱的变化。

要求将xe（n）分成8段，分别采用重叠相加法和重叠保留法。

实验程序：

clc;

n1=0:

1:

3;n2=4:

7;n=0:

518;k=1:

7;

xe=rand（1,512）;

xc1=n1;xc2=8-n2;xc=[xc1,xc2];

xe1=k-k;xe2=[xe1,xe];

%重叠相加法

yn=zeros（1,519）;

yn3=zeros（1,519）;

forj=0:

7

%%重叠相加法%重叠保留法

xj=xe（64*j+1:

64*（j+1））;xj1=xe2（64*j+1:

64*j+71）;

xak=fft（xj,71）;xak1=fft（xj1）;

xck=fft（xc,71）;xck1=fft（xc,71）;

yn1=ifft（xak.*xck）;yn2=ifft（xak1.*xck1）;%每段卷积的结果

temp=zeros（1,519）;temp1=zeros（1,519）;

temp（64*j+1:

64*j+71）=yn1;temp1（64*j+1:

64*j+64）=yn2（8:

71）;

yn=yn+temp;yn3=yn3+temp1;%将每段的卷积结果加到yn中

end;

figure

（1）;

subplot（2,1,1）;plot（n,yn）;

xlabel（'n'）;ylabel（'y（n）'）;title（'xc（n）与xe（n）的线性卷积（重叠相加法）'）;

subplot（2,1,2）;plot（n,abs（fft（yn）））;

xlabel（'k'）;ylabel（'Y（k）'）;axis（[0,600,0,300]）;

figure

（2）;

subplot（2,1,1）;plot（n,yn3）;

xlabel（'n'）;ylabel（'y（n）'）;title（'xc（n）与xe（n）的线性卷积（重叠保留法）'）;

subplot（2,1,2）;plot（n,abs（fft（yn3）））;

xlabel（'k'）;ylabel（'Y（k）的幅频特性'）;axis（[0,600,0,300]）;

运行结果：

重叠相加法：

重叠保留法：

（7）对

做16点DCT变换，取前m个系数

做16点IDCT，得到不同m下的

，求归一化的均方误差

，绘出

随m变化的关系并观察不同的m下的

。

实验结果：

实验三IIR数字滤波器的设计

一、实验目的

（1）掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。

（2）观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲

响应不变法的特点。

（3）熟悉巴特沃思滤波器、