数字信号处理实验报告材料.docx
《数字信号处理实验报告材料.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验报告材料.docx（29页珍藏版）》请在冰豆网上搜索。
数字信号处理实验报告材料
数字信号处理实验讲义
实验一MATLAB简介
实验目的
1．熟悉MATLAB软件的使用方法；
2．MATLAB的绘图功能；
3．用MATLAB语句实现信号的描述及变换。
实验原理
1．在MATLAB下编辑和运行程序
在MATLAB中，对于简单问题可以在命令窗（commandwindows）直接输入命令，得到结果；对于比较复杂的问题则可以将多个命令放在一个脚本文件中，这个脚本文件是以m为扩展名的，所以称之为M文件。
用M文件进行程序的编辑和运行步骤如下：
（1）打开MATLAB，进入其基本界面；
（2）在菜单栏的File项中选择新建一个M文件；
（3）在M文件编辑窗口编写程序；
（4）完成之后，可以在编辑窗口利用Debug工具调试运行程序，在命令窗口查看输出结果；也可以将此文件保存在某个目录中，在MATLAB的基本窗口中的File项中选择RunTheScript，然后选择你所要运行的脚本文件及其路径，即可得出结果；也可以将此文件保存在当前目录中，在MATLAB命令窗口，“>>”提示符后直接输入文件名。
2．MATLAB的绘图功能
plot（x,y）基本绘图函数，绘制x和y之间的坐标图。
figure（n）开设一个图形窗口n
subplot（m,n,N）分割图形窗口的MATLAB函数，用于在一个窗口中显示多个图形，将图形窗口分为m行n列，在第N个窗口内绘制图形。
axis（[a0,b0,a1,b1]）调整坐标轴状态
title（‘’）给图形加题注
xlabel（‘‘）给x轴加标注
ylabel（‘‘）给y轴加标注
grid给图形加网格线
3．信号描述及变换
信号描述及变换包括连续时间信号和离散时间信号内容，详细内容请见课本第1章、第2章。
实验内容
1．试用MATLAB绘制出下列信号的波形：
（Signal1.6）
（1）
；
（2）
（3）
；
（4）
；
（5）
【程序代码】
clearall;closeall;clc;
symst;
x1=exp（-1.5*t）
x2=3*sin（0.5*pi*t）
x3=0.5+0.5*sym（（'sign（t）'））
x4=sym（'heaviside（t）'）+sym（'heaviside（t-1）'）-sym（'2*heaviside（t-2）'）
x5=.5*t*（sym（'heaviside（t）'）-sym（'heaviside（t-4）'））
subplot（2,3,1）;
ezplot（x1）;
axis（[-63-5007000]）;
title（'x1（t）=exp（-1.5t）'）;
gridon
subplot（2,3,2）;
ezplot（x2）;
title（'x2（t）=3sin（0.5¦Ðt）'）;
gridon
subplot（2,3,3）;
fplot（'sign（t）/2+1/2',[-1010],1e-8）;
ezplot（x3,[-1010]）;
axis（[-1010-.21.2]）;
xlabel（'t'）;title（'x3（t）=0.5+0.5sgn（t）'）;
gridon
subplot（2,3,4）;
ezplot（x4,[-13]）;
title（'x4（t）=u（t）+u（t-1）-2u（t-2）'）;
gridon
subplot（2,3,5）;
ezplot（x5,[-26]）;
title（'x5（t）=0.5t[u（t）-u（t-4）]'）;
gridon
subplot（2,3,6）;
axisoff
2．已知连续时间信号（Signal1.7）
，
，
试用MATLAB绘制出下列信号的波形：
（1）
；
（2）
；
（3）
；
（4）
；
（5）
。
【程序代码】
clearall;closeall;clc;
figure
（2）
symst;
x1=（4-t）*（sym（'heaviside（t）'）-sym（'heaviside（t-4）'））
x2=exp（-2*t）*sym（'heaviside（t）'）
x3=sin（2*pi*t）
x4=subs（x1,t,t/2）
x5=subs（x4,t,t-2）
x6=subs（x2,t,-t）
x7=x2+x6
x8=x7*x3
subplot（2,3,1）
text（0,0.9,'x1（t）=（4-t）[u（t）-u（t-4）]'）;
text（0,0.7,'x2（t）=exp（-2t）u（t）'）;
text（0,0.5,'x3（t）=sin（2¦Ðt）'）;
axisoff;boxoff;
subplot（2,3,2）;
ezplot（x5,[-3,10]）;grid
title（'x4（t）=x1（t/2）'）;
subplot（2,3,3）;
ezplot（x5,[-1,12]）;grid
title（'x5（t）=x4（t-2）'）;
subplot（2,3,4）;
ezplot（x6,[-10,5]）;
axis（[-55-0.11.1]）;
grid
title（'x6（t）=x2（-t）'）;
subplot（2,3,5）;
ezplot（x7,[-55]）;
axis（[-55-0.11.1]）;
grid;
title（'x7（t）=x2（t）+x6（t）'）;
subplot（2,3,6）;
tv8=-2.5:
0.05:
2.5;
xv8=subs（x8,tv8）;
plot（tv8,xv8）
axis（[-2.52.5-11]）;
grid
xlabel（'t'）;title（'x8（t）=x7（t）x3（t）'）;
cleartv8xv8
3．列出单位冲激信号、单位阶跃信号、正弦信号的MATLAB表达式，并绘出信号波形。
【程序代码】
clearall;closeall;clc
symst;
x1=sym（'dirac（t）'）;
x2=sym（'Heaviside（t）'）;
x3=sin（t）;
tn=[-6.3:
0.1:
6.3];
xn1=subs（x1,t,tn）;
xn2=subs（x2,t,tn）;
xn3=subs（x3,t,tn）;
plot（tn,xn1,'k',tn,xn2,'r',tn,xn3,'m'）;
grid
xlabel（'t'）;ylabel（'x（t）'）;
legend（'dirac','Heaviside','sin'）
holdon
plot（tn,xn1,'k.',tn,xn2,'r.',tn,xn3,'m.'）
实验二用FFT实现信号的谱分析
实验目的
1．了解FFT在信号谱分析中的作用；
2．了解谱分析的一般步骤和方法。
实验原理
关于信号谱分析的步骤和方法参见教材第3章相关内容。
为了解信号的特点，了解信号频谱分布情况，应该对信号进行谱分析，计算出信号的幅度谱、相位谱和功率谱。
信号的谱分析可以用FFT实现，讨论如下：
1．谱分析中的参数选择；
A若已知信号的最高频率
，为防止混叠，选定采样频率
：
（1）
B根据实际需要，选定频率分辨
，一但选定后，即可确定FFT所需的点数N
（2）
我们希望
越小越好，但
越小，N越大，计算量、存储量也随之增大。
一般取N为2的整次幂，以便用FFT计算，若已给定N，可用补零方法便N为2的整次幂。
C
和N确定后，即可确定所需相应模拟信号
的长度
（3）
分辨率
反比于T，而不是N，在给定的T的情况下，靠减小
来增加N是不能提高分辨率的，因为
为常数
2．谱分析步骤；
A数据准备
（4）
B使用FFT计算信号的频谱
（5）
（6）
C由频谱计算幅度谱
、相位谱
和功率谱
（7）
（8）
（9）
3．实验中用到的一些基本函数简介
y=fft（x,n）；计算n点的FFT。
abs（x）；取绝对值。
angle（z）；取相角。
[Pxx,f]=periodogram（xn,nfft,fs,window）；%周期图谱估计
[Pxx,f]=pwelch（xn,nfft,fs,window,noverlap）；%平均周期图法
Pxx=psd（xn）；功率谱密度
实验内容
1.已知序列x（n）=2sin（0.48πn）+cos（0.52πn）0≤n<100,试绘制x（n）及它的频谱图。
若x（n）=sin（0.56πn）+2cos（0.25πn）,结果又如何？
【程序代码】
clearall;closeall;clc
N=100;
n=0:
N-1;
xn=2*sin（0.48*pi*n）+cos（0.52*pi*n）;
XK=fft（xn,N）;
magXK=abs（XK）;
phaXK=angle（XK）;
subplot（1,2,1）
plot（n,xn）
xlabel（'n'）;ylabel（'x（n）'）;
title（'x（n）N=100'）;
subplot（1,2,2）
k=0:
length（magXK）-1;
stem（k,magXK,'.'）;
xlabel（'k'）;ylabel（'|X（K）|'）;
title（'X（K）N=100'）;
2.对下面信号进行频谱分析，求幅度谱
和相位谱
。
（1）
，
，
，
（2）
，
，
【程序代码】
clearall;closeall;clc
fs1=5000;
dt1=1/fs1;
N1=0.004/dt1;
n1=0:
N1-1;
xn1=0.8.^（n1*dt1）;
Xk1=fft（xn1,N1）;
mag1=2*abs（Xk1）/N1;
pha1=angle（Xk1）;
f1=n1*fs1/N1;
subplot（221）
stem（f1,mag1,'fill'）;
title（'Magnitude1'）
xlabel（'f'）;ylabel（'|Xk
（1）|'）;
subplot（222）
stem（f1,pha1,'fill'）;
title（'Phase1'）
xlabel（'f'）;ylabel（'|Phk
（1）|'）;
fs2=128;
dt2=1/fs2;
N2=18;
n2=0:
N2-1;
xn2=sin（n2*dt2）./（n2*dt2+eps）;
Xk2=fft（xn2,N2）;
mag2=2*abs（Xk2）/N2;
pha2=angle（Xk2）;
f2=n2*fs2/N2;
subplot（223）
stem（f2,mag2,'fill'）;
title（'Magnitude2'）
xlabel（'f'）;ylabel（'|Xk
（2）|'）;
subplot（224）
stem（f2,pha2,'fill'）;
title（'Phase2'）
xlabel（'f'）;ylabel（'|Phk
（2）|'）
3.给定信号
，
，
，现在对
采样，采样点数
，采样频率
=50Hz，设采样序列为
，编写程序计算
的频谱，并绘图；改变采样频率，得到序列
，计算
的频谱，并绘图；增大采样点数，得到序列
，计算
的频谱，并绘图；采样点数N=64，采样频率
=300Hz，在采样点后补零得到新序列
，计算
的频谱，并绘图。
（Signal3.18）
【程序代码】
clearall;closeall;clc
N=16;
n=0:
N-1;
f1=15;f2=18;
fs=50;
dt=1/fs;
xn=sin（2*pi*f1*n*dt）+2*sin（2*pi*f2*n*dt）;
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
figure
（1）;
subplot（121）
plot（f,mag）;grid
title（'Magnitude'）;xlabel（'f'）;ylabel（'|Xk|'）;
subplot（122）
plot（f,pha）;grid
title（'Phase'）;xlabel（'f'）;ylabel（'?
?
?
¡§Xk?
?
'）;
N=16;
n=0:
N-1;
f1=15;f2=18;
fs=36;
dt=1/fs;
xn=sin（2*pi*f1*n*dt）+2*sin（2*pi*f2*n*dt）;
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
figure
（2）;
subplot（221）
plot（f,mag,'r'）;grid
title（'Magnitude（36HzSamplingrequency）'）;
xlabel（'f'）;ylabel（'|Xk|'）;
subplot（222）
plot（f,pha,'r'）;grid
title（'Phase（36HzSamplingrequency）'）;
xlabel（'f'）;ylabel（'?
?
?
¡§Xk?
?
'）;
fs=100;
dt=1/fs;
xn=sin（2*pi*f1*n*dt）+2*sin（2*pi*f2*n*dt）;
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
subplot（223）
plot（f,mag,'r'）;grid
title（'Magnitude（100HzSamplingFrequency）'）;
xlabel（'f'）;ylabel（'|Xk|'）;
subplot（224）
plot（f,pha,'r'）;grid
title（'Phase（100HzSamplingFrequency）'）;
xlabel（'f'）;ylabel（'?
?
?
¡§Xk?
?
'）;
N=64;
n=0:
N-1;
f1=15;f2=18;
fs=50;
dt=1/fs;
xn=sin（2*pi*f1*n*dt）+2*sin（2*pi*f2*n*dt）;
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
figure（3）;
subplot（221）
plot（f,mag,'m'）;grid
title（'Magnitude（64SamplingPoints）'）;xlabel（'f'）;ylabel（'|Xk|'）;
subplot（222）
plot（f,pha,'m'）;grid
title（'Phase（64SamplingPoints）'）;xlabel（'f'）;ylabel（'?
?
?
¡§Xk?
?
'）;
N=512;
n=0:
N-1;
dt=1/fs;
xn=sin（2*pi*f1*n*dt）+2*sin（2*pi*f2*n*dt）;
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
subplot（223）
plot（f,mag,'m'）;grid
title（'Magnitude（512SamplingPoints）'）;xlabel（'f'）;ylabel（'|Xk|'）;
subplot（224）
plot（f,pha,'m'）;grid
title（'Phase（512SamplingPoints）'）;xlabel（'f'）;ylabel（'?
?
?
¡§Xk?
?
'）;
N=64;
n=0:
N-1;
f1=15;f2=18;
fs=300;
dt=1/fs;
xn0=sin（2*pi*f1*n*dt）+2*sin（2*pi*f2*n*dt）;
xn=[xn0,zeros（1,64）];
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
figure（4）;
subplot（221）
plot（f,mag,'g'）;grid
title（'Magnitude（filledwith64zeros）'）;xlabel（'f'）;ylabel（'|Xk|'）;
subplot（222）
plot（f,pha,'g'）;grid
title（'Phase（filledwith64zeros）'）;xlabel（'f'）;ylabel（'?
?
?
¡§Xk?
?
'）;
xn=[xn0,zeros（1,448）];
Xk=fft（xn,N）;
mag=2*abs（Xk）/N;
pha=angle（Xk）;
f=fs/N*n;
subplot（223）
plot（f,mag,'g'）;grid
title（'Magnitude（filledwith448zeros）'）;xlabel（'f'）;ylabel（'|Xk|'）;
4.试求下列差分方程所描述的输出序列
的功率谱并作图。
（a）
（b）
（c）
式中，
是方差为
（例如，
=1/12）的白噪声。
【程序代码】
clearall;closeall;clc
nfft=512;
window=hanning（256）;
noverlap=128;
w=randn（1,1000）;
fori=1:
1000
ifi<=2
x（i）=randn（1,1）;
else
x（i）=-0.81*x（i-2）+w（i）-w（i-1）;
end
end
[Pxxf]=psd（x,nfft,2,window,noverlap,'mean'）;
figure
（1）
plot（f,10*log10（Pxx））;grid
xlabel（'f'）;ylabel（'PowerSpectrum,dB'）;
title（'（a）:
x（n）=-0.81x（n-2）+w（n）-w（n-1）'）
w=randn（1,1000）;
fori=1:
1000
ifi<=2
x（i）=randn（1,1）;
else
x（i）=w（i）-w（i-2）;
end
end
[Pxxf]=psd（x,nfft,2,window,noverlap,'mean'）;
figure
（2）
plot（f,10*log10（Pxx）,'m'）;grid
xlabel（'f'）;ylabel（'PowerSpectrum,dB'）;
title（'（b）:
x（n）=w（n）-w（n-2）'）
w=randn（1,1000）;
fori=1:
1000
ifi<=2
x（i）=randn（1,1）;
else
x（i）=-0.81*x（i-2）+w
（2）;
end
end
[Pxxf]=psd（x,nfft,2,window,noverlap）;
figure（3）
plot（f,10*log10（Pxx）,'r'）;grid
xlabel（'f'）;ylabel（'PowerSpectrum,dB'）;
title（'（c）:
x（n）=-0.81x（n-2）+w（n）'）
5.一序列
是由两个频率相距为
的模拟信号采样得来的，即
n=0,1,…,15
已知序列长度N=16，试采用周期图法，应用DFT分别计算当
=0.06及
=0.01时的功率谱估计，并通过作图说明从功率谱估计的分布是否能分辨出这两个正弦信号的真实频谱？
若N=64又有什么变化？
【程序代码】
clearall;closeall;clc
nfft=512;window=hanning（256）;noverlap=128;df=[0.060.01];
N=16;n=0:
N-1;
x1=sin（2*pi*0.135*n）+cos（2*pi*（0.135+df
（1））*n）;
x2=sin（2*pi*0.135*n）+cos（2*pi*（0.135+df
（2））*n）;
[Pxx1f1]=psd（x1,nfft,1,window,noverlap,'mean'）;
[Pxx2f2]=psd（x2,nfft,1,window,noverlap,'mean'）;
figure
（1）
plot（f1,Pxx1）;grid
xlabel（'f'）;ylabel（'PowerSpectrum'）;
title（'sin2?
?
（0.135）n+cos2?
?
（0.135+0.06）nN=16'）
figure
（2）
plot（f2,Pxx2）;grid
xlabel（'f'）;ylabel（'PowerSpectrum'）;
title（'（sin2?
?
（0.135）n+cos2?
?
（0.135+0.01）nN=16'）
N=64;n=0:
N-1;
x1=sin（2*pi*0.135*n）+cos（2*pi*（0.135+df
（1））*n）;
x2=sin（2*pi*0.135*n）+cos（2*pi*（0.135+df
（2））*n）;
[Pxx1f1]=psd（x1,nfft,1,window,noverlap,'mean'）;
[Pxx2f2]=psd（x2,nfft,1,window,noverlap,'mean'）;
figure（3）
plot（f1,Pxx1,'m'）;grid
xlabel（'f'）;ylabel（'PowerSpectrum'）;
title（'sin2?
?
（0.135）n+cos2?
?
（0.135+0.06）nN=64'）
figure（4）;plot（f2,Pxx2,'m'）;grid
xlabel（'f'）;ylabel（'PowerSpectrum'）;
title（'（sin2?
?
（0.135）n+cos2?
?
（0.135+0.01）nN=64'）
6.用MATLAB产生256点白噪声序列，应用Welch法估计其功率谱，每段长64点，重叠32点，输出平均后的功率谱图以及对256点一次求周期图的功率谱图。
【程序代码】
N=256;
x=rand（1,N）;
window=hanning（64）;
noverlap=32;nfft=128;
[Pxx1f1]=periodogram（x）;
figure
（1）
plot（f1/pi,10*log10（Pxx1））;grid
xlabel（'f'）;ylabel（'PowerSpectrum,dB'）;
title（'¶Ô256µãÒ»´ÎÇóÖÜÆÚͼµÄ¹¦ÂÊÆ×ͼ'）
[Pxx2f2]=psd（x,nfft,2,window,noverlap,'line'）;
figure
（2）
plot（f2,10*log10（Pxx2）,'
升级会员 