数字信号处理MATLAB实验

1、适用于本科生所有的实验报告（印制实验报告册除外）；

2、专业填写为专业全称，有专业方向的用小括号标明；

3、格式要求：

①用A4纸双面打印（封面双面打印）或在A4大小纸上用蓝黑色水笔书

写。

②打印排版：正文用宋体小四号，1.5倍行距，页边距采取默认形式（上

下2.54cm，左右2.54cm，页眉1.5cm，页脚1.75cm）。字符间距为默认

值（缩放100%，间距：标准）；页码用小五号字底端居中。

③具体要求：

题目（二号黑体居中）；

摘要（“摘要”二字用小二号黑体居中，隔行书写摘要的文字部分，小4号宋体）；

关键词（隔行顶格书写“关键词”三字，提炼3-5个关键词，用分号隔

开，小4号黑体)；

正文部分采用三级标题；

第1章××(小二号黑体居中，段前0.5行)

1.1 ×××××小三号黑体×××××（段前、段后0.5行）

1.1.1小四号黑体（段前、段后0.5行）

参考文献（黑体小二号居中，段前0.5行），参考文献用五号宋体，参照《参考文献著录规则（GB/T 7714－2005）》。

实验一 利用MATLAB 求解白噪声、

相关函数

实验内容及要求：

一、MATLAB 文件学习（1学时）

学习《1.7节 与本章内容相关的MATLAB 文件》，包括（1）rand.m 、（2）randn.m 、（3）sinc.m 、（4）xcorr.m ，并熟练掌握其功能、调用格式及其参数的含义。

二、上机练习（3学时）

练习第一章课后习题的*1.14、*1.15、*1.16，并将程序代码及相关结果（包括文字、图形）写入实验报告。

题目一：

习题1.14（1）：实验内容：把连续时间信号()x t 转换成离散时间信号()s x nT ，并绘出()s x nT 的图形。1/s s f T =为抽样频率。

（1）工频信号：10()sin(2)x t A f t π=，220A =，050f Hz =

程序代码：

clear

t=-0.05:0.001:0.05; A=220;

fo=50;fs=5;

x=A*sin(2*pi*fo*t);

Ts=1/fs;n=-0.05:0.001:0.05; xl=A*sin(2*pi*fo*n*Ts); subplot(211)

stem(n,xl,'.');grid on;xlabel('n');ylabel('xl'); subplot(212)

plot(t,x);grid on;xlabel('t');ylabel('xl');

图1.1 fs=5Hz时()

x nT的图形

s

图1.2 fs=10Hz时()

x nT的图形

s

图1.3 fs=25Hz时()

x nT的图形

s

图1.4 fs=50Hz时()

x nT的图形

s

图1.5 fs=100Hz时()

x nT的图形

s

图1.6 fs=250Hz时()

x nT的图形

s

题目二：

习题1.14（3）：实验内容：把连续时间信号()x t 转换成离散时间信号()s x nT ，并绘出()s x nT 的图形。1/s s f T =为抽样频率。

（3）谐波信号：x 3(t )=∑A i sin (2πf 0it )3i=1，其中A 1=1，A 2=0.5，A 3=0.2，f 0=5Hz

程序代码：

clear

t=-0.05:0.001:0.05;

A1=1;A2=0.5;A3=0.2;fo=5;fs=25;

xl=A1*sin(2*pi*fo*t)+A2*sin(2*pi*fo*2*t)+A3*sin(2*pi*fo*3*t);

Ts=1/fs;n=-0.05:0.001:0.05;

x=A1*sin(2*pi*fo*n*Ts)+A2*sin(2*pi*fo*2*n*Ts)+A3*sin(2*pi*fo*3*n*Ts); subplot(211)

stem(n,x'.');grid on;xlabel('n');ylabel('x'); subplot(212)

plot(t,xl);grid on;xlabel('t');ylabel('xl');

实验结果：

图2.1 fs=2.5Hz时()

x nT的图形

s

图2.2 fs=5Hz时()

x nT的图形

s

图2.3 fs=10Hz时()

x nT的图形

s

图2.4 fs=25Hz时()

x nT的图形

s

题目三：

习题1.14（5）：实验内容：把连续时间信号()x t 转换成离散时间信号()s x nT ，并绘出()s x nT 的图形。1/s s f T =为抽样频率。

（5）sinc 函数：x 3(t )=sin (Ωt )

Ωt ⁄,其中Ω=2πf ，f =10Hz

程序代码：

（文字、图形等）

clear

f=10;omega=2*pi*f; t=-0.5:0.001:0.5; x=omega*t; y=sinc(x); subplot(212)

plot(t,y);grid on ; fs=5;

n=-0.5:0.01:0.5; xl=omega*n*(1/fs); yl=sinc(xl); subplot(211)

stem(n,yl);grid on ;

实验结果：

图3.1 fs=5Hz时()

x nT的图形

s

图3.2 fs=10Hz时()

x nT的图形

s

图3.3 fs=20Hz时()

x nT的图形

s

图3.4 fs=50Hz时()

x nT的图形

s题目四：

习题1.15:实验内容：令x(n)=Asin(ωn)+u(n),其中ω=π/16，u(n)是白噪声。（1）使用MATLAB中的有关文件，产生均值为0功率P=0.1的均匀分布白噪声u(n)，画出其图形，并求u(n)的自相关r

(m),画出r u(m)的波形。

u

（2）欲使x(n)的信噪比为10dB，试决定A的数值，并画出x(n)的图形及其自相关函数r

(m)的图形。

x

程序代码：

p=0.1;N=50000;n=(1:100);

u=rand(1,N);u=u-mean(u);

a=sqrt(12*p);ul=u*a;

power_ul=var(ul)/N;

subplot(221)

plot(ul(1:100));grid on;xlabel('u');ylabel('ul');

title('白噪声u(n)的图形');

ru=xcorr(ul(1:100));

subplot(222)

plot(ru);grid on;xlabel('ul');ylabel('ru');

title('u(n)的自相关函数ru (m)的波形');

Ps=1;Pn=0.01;w=pi/16;

snr=10*log10(Ps/Pn);

A=sqrt(Ps*2);

x=A*sin(w*n)+ul(1:100);

rx=xcorr(x(1:100));

subplot(223)

plot(n,x);grid on;xlabel('n');ylabel('x');

title('x(n)的图形');

subplot(224)

plot(rx);grid on;xlabel('x');ylabel('rx');

title('x(n)的自相关函数rx (m)的图形');

实验结果：

题目五：

习题1.16:实验内容：表中给出的是从1770——1869年这100年间每一年12个月所记录到的太阳黑子出现的次数的平均值。

（1）输出该数据的图形。

（2）对该数据做自相关，输出其自相关函数的图形，观察太阳黑子活动的周期（取M=32）。

（3）将该数据除去均值。再重复（2）的内容，比较除去均值前、后对做自相关的影响。

程序代码：

clear

A=[101,82,66,35,31,7,20,92,154,125,85,68,38,23,10,24,83,132,131,118,90,67,60,47,41,21,16,6,4,7,14,34,45,43 ,48,42,28,10,8,2,0,1,5,12,14,35,46,41,30,24,16,16,4,2 ,8,17,36,50,62,67,71,48,28,8,13,57,122,138,103,86,63, 37,24,11,15,40,62,98,124,96,66,,54,39,21,7,4,23,55, 94,96,77,59,44,4,30,16,7,37,74];

B=1770:1869;

subplot(311)

stem(B,A,'.');grid on;xlabel('时间');

ylabel('太阳黑子出现次数');

M=32;

ra=xcorr(A,M);

subplot(312)

plot(ra);grid on;ylabel('自相关系数');

title('除去均值前的自相关函数图形');

A=A-mean(A);

rb=xcorr(A,M);

subplot(313)

plot(rb);grid on;ylabel('自相关系数');

title('除去均值后的自相关函数图形');

实验结果：

实验二利用MATLAB文件求解离散时间系

统的输出、单位抽样响应、

频率响应、极零图

实验内容及要求：

一、MATLAB文件学习（1学时）

学习《2.11节与本章内容相关的MATLAB文件》，包括（1）conv.m、（2）filter.m、（3）impz.m、（4）frqz.m、(5)zplane.m，并熟练掌握其功能、调用格式及其参数的含义。

二、上机练习（3学时）

练习第一章课后习题的*2.17，并将程序代码及相关结果（包括文字、图形）写入实验报告。

题目一：

习题2.17：实验内容：给定系统H(z)=−0.2z/(z2+0.8)。

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

（2）求出并绘出该系统的单位抽样响应ℎ(n)。

（3）令x(n)=u(n)，求出并绘出系统的单位阶跃响应y(n)。

程序代码：

clear

x=ones(100);

t=1:100;

b=[0,-0.2,0];

a=[1,0,0.8];

[h,t]=impz(b,a,40);

y=filter(b,a,x);

[H,w]=freqz(b,a,256,'whole',1);Hr=abs(H);

Hphase=angle(H);

Hphase=unwrap(Hphase);

subplot(221)

stem(t,h);title('单位抽样响应');

xlabel('t');ylabel('h');grid on; subplot(222)

plot(y);title('阶跃响应');

xlabel('x');ylabel('y');grid on; subplot(223)

plot(w,Hr);title('幅频响应');

xlabel('w');ylabel('Hr');grid on; subplot(224)

plot(w,Hphase);title('相频响应');

xlabel('w');ylabel('Hphase');grid on; 实验结果：

实验三：利用MATLAB求解系统滤波

实验内容及要求：

一、MATLAB文件学习（1学时）

学习《第三章3.9节与本章内容相关的MATLAB文件》，熟练掌握其功能、调用格式及其参数的含义。

二、上机练习（3学时）

练习第三章课后习题的3.20 *3.21，并将程序代码及相关结果（包括文字、图形）写入实验报告。

题目一：

3.20 进一步研究参数T和τ对傅里叶图形的影响。例如，分别令τ=0.2T,0.1T及0.05T，试用MATLAB求出并画出傅里叶系数图，并分析这些参数变化对图形的规律。

程序代码：

clear

A=5;T=1;

tao=0.2*T;

k=-30:30;

x=k*pi*tao;

y=sinc(x);

z=A*tao/T*y;

subplot(311);stem(k,z,'.');title('τ=0.1T时');

xlabel('k');ylabel('z');grid on

tao=0.1*T;

k=-30:30;

x=k*pi*tao;

y=sinc(x);

z=A*tao/T*y;subplot(312);stem(k,z,'.');title('τ=0.2T时');

xlabel('k');ylabel('z');grid on

tao=0.05*T;

k=-30:30;

x=k*pi*tao;

y=sinc(x);

z=A*tao/T*y;

subplot(313);stem(k,z,'.');title('τ=0.05T时'); xlabel('k');ylabel('z');grid on

实验结果：

题目二：

3.21 设有一长序列x(n)

x(n){n

5⁄ 0≤n≤50

20−n5⁄ 500 其他

令x(n)通过一离散系统，其单位抽样响应

ℎ(n){1

2n

⁄ 0≤n≤2 0 其他

试编一主程序用叠加相加法实现该系统对x(n)的滤波，并画出输出y(n)的图形。程序代码：

clear

n=0:1:100;

x=zeros(1,length(n));

h=zeros(1,length(n));

y=zeros(1,length(n));

for m=1:length(n)

if n(m)>=0&&n(m)<=50

x(m)=n(m)/5;

elseif n(m)>50&&n(m)<=99

x(m)=20-n(m)/5;

else

x(m)=0;

end

end

for m=1:length(n)

if n(m)>=0&&n(m)<=2

h(m)=0.5^n(m);

else

h(m)=0;

end

end

y=fftfilt(h,x);

plot(y);grid on

xlabel('n')

ylabel('y(n)')

title('重叠相加法');

实验结果：

实验四快速傅里叶变换

实验内容及要求：

一、MATLAB文件学习（2学时）

学习《4.7节与本章内容相关的MATLAB文件》，包括fft.m，czt.m并熟练掌握其功能、调用格式及其参数的含义。

二、上机练习（2学时）

练习第三章课后习题的*4.10并将程序代码及相关结果（包括文字、图形）写入实验报告。

题目一：

4.10 CZT算法的研究。给定信号

3

x(t)=∑sin (2πf i t)

i=1

已知f1=10.8Hz,f2=11.75Hz,f3=12.55Hz，令,f4=40Hz，对x(t)抽样后得x(n)，又令N=。

=（1）调用MATLAB中的czt.m，可求出X(k)及幅度谱，此时△f=f s

N

0.625Hz,小于(f2−f1)及(f3−f2)，观察三个谱峰的分辨情况。

（2）在x(n)后分别补3N个零、7N个零、15N个零，在做DFT，观察补零的效果。

（3）调用MATLAB中的文件czt.m，按如下两组参数赋值。

参数1：f s=40Hz，N=，M=50，f0=9Hz，△f=0.2Hz；

参数2：f s=40Hz，N=，M=60，f0=8Hz，△f=0.12Hz。

分别求X(k)，k=0，1，⋯，M−1，画出其幅度谱，并和（1）（2）的结果相比较。

程序代码：

（1）

clear

N=;f1=10.8;f2=11.75;f3=12.55;fs=40;A=1;

n=0:N-1;

t=2*pi*n/fs;

x=sin(f1*t)+sin(f2*t)+sin(f3*t);

M=N;W=exp(-1i*2*pi/N);

stepf=fs/N;

n1=0:stepf:fs/2-stepf;

Y=czt(x,M,W,A);

subplot(2,2,k); plot(n1, abs(Y(1: N/2))); grid on; xlabel('f/Hz');ylabel('X(k)');

（2）

clear

N=;f1=10.8;f2=11.75;f3=12.55;fs=40;A=1;

n=0:N-1;

t=2*pi*n/fs;

x=sin(f1*t)+sin(f2*t)+sin(f3*t);

M=N;W=exp(-1i*2*pi/N);

stepf=fs/N;

n1=0:stepf:fs/2-stepf;

Y2=abs(fft(x));

subplot(2,2,1); plot(n1, abs(Y2(1: N/2))); grid on;

xlabel('f/Hz');ylabel('X(k)的幅度谱');title('用FFT计算DFT')

x(N:4*N-1)=0; stepf=fs/(4*N);

n1=0:stepf:fs/2-stepf;

Y2=abs(fft(x));

subplot(2,2,2);plot(n1,abs(Y2(1:2*N)));grid on;

xlabel('f/Hz');ylabel('X(k)的幅度谱');title('补3N个零');

x(N:8*N-1)=0; stepf=fs/(8*N);

n1=0:stepf:fs/2-stepf;

Y2=abs(fft(x));

subplot(2,2,3);plot(n1,abs(Y2(1:4*N)));grid on;

xlabel('f/Hz');ylabel('X(k)的幅度谱');title('补7N个零');x(N:16*N-1)=0; stepf=fs/(16*N);

n1=0:stepf:fs/2-stepf;

Y2=abs(fft(x));

subplot(2,2,4);plot(n1,abs(Y2(1:8*N)));grid on;

xlabel('f/Hz');ylabel('X(k)的幅度谱');title('补15N个零'); （3）

clear

N=;f1=10.8;f2=11.75;f3=12.55;fs=40;A=1;n=0:N-1;

t=2*pi*n/fs;

x=sin(f1*t)+sin(f2*t)+sin(f3*t);

N1=[1,3,7,15];

fs=40;N=;M=50;f0=9;stepf=0.2;A=exp(1i*2*pi*f0/fs);

%参数一

for k=1:4

N=N1(k)*;W=exp(-1i*2*stepf/fs);

n1=f0:stepf:f0+(M-1)*stepf;

Y=czt(x,M,W,A);

subplot(2,2,k); plot(n1, abs(Y)); grid on;

xlabel('f/Hz');ylabel('X1(e^j^w)');

end

fs=40;N=;M=60;f0=8;stepf=0.12;A=exp(1i*2*pi*f0/fs); %参数二

for k=1:4

N=N1(k)*;W=exp(-1i*2*stepf/fs);

n1=f0:stepf:f0+(M-1)*stepf;

Y=czt(x,M,W,A);

subplot(2,2,k); plot(n1, abs(Y)); grid on;

xlabel('f/Hz');ylabel('X1(e^j^w)');

end

实验结果：

图一：幅度谱与补零

图二：参数一与补零

图三：参数三与补零