数字信号管理目标实验报告.docx
《数字信号管理目标实验报告.docx》由会员分享，可在线阅读，更多相关《数字信号管理目标实验报告.docx（22页珍藏版）》请在冰点文库上搜索。
//
数字信号处理实验讲义
前言 2
实验一 MATLAB简介 3
实验二 用FFT实现信号的谱分析 5
实验三 IIR数字巴特沃思滤波器的设计 8
实验四 FIR数字滤波器的设计 9
实验一 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
（1） x（t）=e-1.5t；
（2） x2（t）=3sin（0.5pt）
（3） x3（t）=0.5+0.5sgn（t）；
（4）
（5）（5）
x4（t）=u（t）+u（t-1）-2u（t-2）；
x（t）=t[u（t）-u（t-4）]
5 2
【程序代码】
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
3
2.已知连续时间信号（Signal1.7）
x1（t）=（4-t）[u（t）-u（t-4）]，x2
（t）=e-2tu（t），x
（t）=sin（2pt）
试用MATLAB绘制出下列信号的波形：
（1） x4（t）=x1（t/2）；
（2） x5（t）=x4（t-2）；
（3） x6（t）=x2（-t）；
（4） x7（t）=x2（t）+x6（t）；
（5） x8（t）=x7（t）×x3（t）。
【程序代码】
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*x3subplot（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]）;
gridxlabel（'t'）;title（'x8（t）=x7（t）x3（t）'）;cleartv8xv8
3.列出单位冲激信号、单位阶跃信号、正弦信号的MATLAB表达式，并绘出信号波形。
【程序代码】
clearall;closeall;clcsymst;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'）holdonplot（tn,xn1,'k.',tn,xn2,'r.',tn,xn3,'m.'）
实验二 用FFT实现信号的谱分析
实验目的
1.了解FFT在信号谱分析中的作用；
2.了解谱分析的一般步骤和方法。
实验原理
关于信号谱分析的步骤和方法参见教材第3章相关内容。
为了解信号的特点，了解信号频谱分布情况，应该对信号进行谱分析，计算出信号的幅度谱、相位谱和功率谱。
信号的谱分析可以用FFT实现，讨论如下：
1.谱分析中的参数选择；
A 若已知信号的最高频率fc，为防止混叠，选定采样频率fs：
fs³2fc
（1）
B根据实际需要，选定频率分辨Df，一但选定后，即可确定FFT所需的点数N
N=fs/Df
（2）
我们希望Df越小越好，但Df越小，N越大，计算量、存储量也随之增大。
一般取N为2的整次幂，以便用FFT计算，若已给定N，可用补零方法便N为2的整次幂。
C fs和N确定后，即可确定所需相应模拟信号x（t）的长度
T=N/fs=NTs
（3）
分辨率Df反比于T，而不是N，在给定的T的情况下，靠减小Ts来增加N是不能提高分辨率的，因为T=NTs为常数
2.谱分析步骤；A数据准备
x（n）=xa（t）t=nT=xa（nT）
B使用FFT计算信号的频谱
（4）
X（k）=åx（n）W
N-1
kn
N
（5）
n=0
X（k）=Xr（k）+jXi（k）
C由频谱计算幅度谱X（k）、相位谱qk和功率谱G（k）
（6）
X2（k）+X2（k）
r
i
X（k）=
q=arctanXi（k）
（7）
（8）
X
r
k （k）
G（k）=
X（k）2=X2（k）+X2（k）
（9）
r i
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.对下面信号进行频谱分析，求幅度谱X（k）和相位谱q（k）。
1
（1）x（t）=at，a=0.8，0£t£4ms，f
（2）x2（t）=sint/t，T=0.125s，N=16
max
=400Hz
【程序代码】
clearall;closeall;clcfs1=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. 给定信号x（t）=sin（2pf1t）+2sin（2pf2t），f1=15Hz，f1=18Hz，现在对
x（t）采样，采样点数N=16，采样频率fs=50Hz，设采样序列为x（n），编写程序计算
x（n）的频谱，并绘图；改变采样频率，得到序列x1（n），计算x1（n）的频谱，并绘图；增
大采样点数，得到序列x2（n），计算x2（n）的频谱，并绘图；采样点数N=64，采样频率
fs=300Hz，在采样点后补零得到新序列x3（n），计算x3（n）的频谱，并绘图。
（Signal3.18）
【程序代码】
clearall;closeall;clcN=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）;gridtitle（'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.试求下列差分方程所描述的输出序列x（n）的功率谱并作图。
（a）x（n）=-0.81x（n-2）+w（n）-w（n-1）
（b）x（n）=w（n）-w（n-2）
（c）x（n）=-0.81x（n-2）+w（n）
式中，w（n）是方差为s2（例如，s2=1/12）的白噪声。
w w
【程序代码】
clearall;closeall;clcnfft=512;window=hanning（256）;noverlap=128;w=randn（1,1000）;
fori=1:
1000ifi<=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））;gridxlabel（'f'）;ylabel（'PowerSpectrum,dB'）;title（'（a）:
x（n）=-0.81x（n-2）+w（n）-w（n-1）'）w=randn（1,1000）;
fori=1:
1000ifi<=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'）;gridxlabel（'f'）;ylabel（'PowerSpectrum,dB'）;title（'（b）:
x（n）=w（n）-w（n-2）'）w=randn（1,1000）;
fori=1:
1000ifi<=2
x（i）=randn（1,1）;else
x（i）=-0.81*x（i-2）+w
（2）;
endend
[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.一序列x（n）是由两个频率相距为Df的模拟信号采样得来的，即
x（n）=sin2p（0.135）n+cos2p（0.135+Df）n
n=0,1,…,15
已知序列长度N=16，试采用周期图法，应用DFT分别计算当Df=0.06及Df=0.01时的功
率谱估计，并通过作图说明从功率谱估计的分布是否能分辨出这两个正弦信号的真实频谱？
若N=64又有什么变化？
【程序代码】
clearall;closeall;clcnfft=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,'mea
升级会员 