数字信号处理实验指导书.docx
《数字信号处理实验指导书.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书.docx（14页珍藏版）》请在冰豆网上搜索。
数字信号处理实验指导书
数字信号处理
实验一：
信号的时域分析——波形的产生及MATLAB实践---------1
实验二：
系统响应及系统稳定性--------------------------------------------3
实验三：
用FFT对信号做频谱分析----------------------------------------5
实验四：
IIR数字滤波器设计与软件实现---------------------------------8
实验一：
信号的时域分析——波形的产生及MATLAB实现
一、实验目的
1、熟悉单位脉冲序列、单位阶跃序列、矩形序列和指数序列的产生。
2、掌握利用MATLAB画图函数显示信号波形的方法。
3、熟悉离散时间信号的翻转、移位和尺度变换等运算
二、实验内容
1、利用MATLAB分别产生并显示下列序列
（1）单位脉冲序列
（2）单位阶跃序列
（3）矩形序列
2、利用MATLAB分别产生并显示下列周期波形。
（1）幅度为1，脉冲宽度占空比duty=1/2，周期T=6的周期方波脉冲信号；
（2）幅度为±1，周期T=4的三角波脉冲.
三、MATLAB实验
上机实验时，用WINDOWS操作系统的记事本或用MATLAB提供的编辑器，逐行输入下列语句，然后以*.m存盘，其中m为后缀名。
所有MATLAB的程序都必须以m作为后缀名保存。
在MATLAB命令窗口，改变当前路径至存有*.m的目录，可用dir命令检验。
最后，在MATLAB命令窗口敲入文件名并回车，则可得到运行结果。
例一：
单位脉冲序列的产生
n=[-3:
3];xn=[（n-0）==0];序列的起点为-3，终点为3，在n=0处生成一个单位脉冲；
subplot（3,2,1）;stem（n,xn,'.'）;stem：
离散序列绘图指令；
line（[-3,3],[0,0]）创建直线；
axis（[-3,3,0,1]）;控制坐标轴比例和外观；
xlabel（'n'）;ylabel（'x（n）'）X轴和Y轴标注；
例二：
单位阶跃序列的产生
n=[-10:
10];xn=[（n-0）>=0];序列的起点为-10，终点为10，在n=0处生成一个单位阶跃；
subplot（3,2,1）;stem（n,xn,'.'）;line（[-10,10],[0,0]）
axis（[-10,10,0,1]）;
xlabel（'n'）;ylabel（'x（n）'）
例三：
矩形序列的产生
clear,closeall清除所有程序；
n=[-10:
10];xn1=[（n-0）>=0];
xn2=[（n-5）>=0];定义两个阶跃序列；
xn=xn1-xn2;两个阶跃序列之差得到矩形序列；
subplot（3,2,1）;stem（n,xn,'.'）;line（[-10,10],[0,0]）
axis（[-10,10,0,1]）;
xlabel（'n'）;ylabel（'x（n）'）
例四：
幅度为1，脉冲宽度占空比duty=50，周期T=5的周期方波脉冲信号的产生。
T=6,t=-2*T:
0.01:
2*T;duty=50;duty—信号为正的区域在一个周期内所占的百分比；
x=square（t,duty）;
subplot（2,2,1）
plot（t,x）
title（'连续时间函数-周期方波脉冲'）title—标出图名；
axis（[-10,10,-1.2,1.2]）
line（[-10,10],[0,0]）
例五：
幅度为±1，周期T=4的三角波脉冲信号的产生
T=4;t=-2.5*T:
0.01:
2.5*T;
x=sawtooth（t-2*T/3,0.5）；sawtooth（x,width）：
三角波函数，width为该位置的横坐标和周期的比值（0-1的标量）。
width=1时产生锯齿波，width=0.5时产生标准对称三角波；
subplot（3,3,1）;
plot（t,x）;
title（'周期三角波脉冲，T=4'）
axis（[-10,10,-1.2,1.2]）;
line（[-10,10],[0,0]）
实验二：
系统响应及系统稳定性
一、实验目的
（1）掌握求系统响应的方法。
（2）掌握时域离散系统的时域特性。
（3）分析观察及检验系统的稳定性。
二、实验原理与方法（略）
三、实验内容及步骤（略）
四、实验参考程序
1、系统响应及系统稳定性
一个低通滤波器的差分方程为
y（n）=0.05x（n）+0.05x（n-1）+0.9y（n-1）
输入信号为
，
1分别求他们的系统响应
和
，并画出其波形。
2求出系统的单位脉冲响应，画出波形。
clear;closeall;
N=64;n=[0:
N-1];
A=[1,-0.9];B=[0.05,0.05];系统差分方程系数向量；
x=[n==0];产生deta（n）；
x1n=[11111111zeros（1,56）];产生
x2n=ones（1,64）;产生
hn=filter（B,A,x）;filter求系统对序列x的响应；
subplot（2,2,1）;
stem（n,hn,'.'）;
line（[0,64],[0,0]）;
axis（[0,64,0,0.1]）;
title（'（a）系统单位脉冲响应h（n）'）
y1n=filter（B,A,x1n）求系统对序列x1n的响应；
subplot（2,2,2）;
stem（n,y1n,'.'）;line（[0,64],[0,0]）;
axis（[0,64,0,0.6]）;
title（'（b）系统对R8（n）的响应y1（n）'）
y2n=filter（B,A,x2n）;求系统对序列x2n的响应；
subplot（2,2,4）;
stem（n,y2n,'.'）;line（[0,64],[0,0]）;
axis（[0,64,0,1]）;
title（'（c）系统对u（n）的响应y2（n）'）
其中求hn也可用impz来实现，改为hn=impz（B,A,64），就可以不用定义序列x。
2、给定系统的单位脉冲响应为
用线性卷积法求
分别对系统
和
的输出响应
和
，并画出波形。
closeall;clearall
x1n=[11111111];N=8;M1=10;M2=4;
h1n=ones（1,10）;
h2n=[12.52.51];
y1n=conv（x1n,h1n）;conv求卷积；
y2n=conv（x1n,h2n）;
L1=N+M1-1;L2=N+M2-1;
n1=0:
L1-1;n2=0:
L2-1;
subplot（2,2,1）;
stem（n1,y1n,'.'）;
line（[0,20],[0,0]）;
axis（[0,20,0,10]）;
title（'（a）h1（n）与R8（n）的卷积'）;
subplot（2,2,2）;
stem（n2,y2n,'.'）;
line（[0,15],[0,0]）;
axis（[0,15,0,10]）;
title（'（b）h2（n）与R8（n）的卷积'）
思考：
1、如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应？
如何求？
可用分段线性卷积法来求解。
2、如果信号经过低通滤波器，信号的高频分量被滤掉，时域信号会有何变化？
当高频分量被滤掉，时域信号的变化减缓，在有阶跃处附近产生过渡带。
因此，当输入巨星序列时，输出序列的开始和终了都产生了明显的过渡带。
见第一个实验波形。
实验三：
用FFT对信号做频谱分析
一、实验目的（略）
二、实验原理（略）
三、实验内容及步骤（略）
四、相关程序
（1）对以下序列进行谱分析
选择变换区间N=8和N=16进行谱分析，分别画出幅频特性曲线。
clearall;closeall;
x1n=ones（1,4）;
N1=8;N2=16;k1=0:
N1-1;k2=0:
N2-1;w1=k1*2/N1;w2=k2*2/N2;k1和k2是频域取样点的取值区间，w1和w2是将k值转换为归一化频率；
M=8;xa=1:
（M/2）;xb=（M/2）:
-1:
1;x2n=[xa,xb];构造x2n
x3n=[xb,xa];
X1k8=fft（x1n,8）;
X1k16=fft（x1n,16）;
X2k8=fft（x2n,8）;
X2k16=fft（x2n,16）;
X3k8=fft（x3n,8）;
X3k16=fft（x3n,16）;
subplot（2,2,1）;stem（w1,abs（X1k8）,'.'）;
title（'（1a）8点DFT[x1（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k8））]）;
subplot（2,2,2）;stem（w2,abs（X1k16）,'.'）;
title（'（1b）16点DFT[x1（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k16））]）;
figure
（2）创建图形窗
subplot（2,2,1）;stem（w1,abs（X2k8）,'.'）;
title（'（2a）8点DFT[x2（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k8））]）;
subplot（2,2,2）;stem（w2,abs（X2k16）,'.'）;
title（'（2b）16点DFT[x2（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k16））]）;
subplot（2,2,3）;stem（w1,abs（X3k8）,'.'）;
title（'（3a）8点DFT[x3（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k8））]）;
subplot（2,2,4）;stem（w2,abs（X3k16）,'.'）;
title（'（3b）16点DFT[x3（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X3k16））]）;
（2）对以下周期序列进行谱分析
选择FFT的变换区间N为8和16两种情况进行谱分析。
分别画出幅频特性曲线。
clearall;closeall;
N1=8;N2=16;k1=0:
N1-1;k2=0:
N2-1;w1=k1*2/N1;w2=k2*2/N2;n=0:
15;
x1n=cos（pi*n/4）;
x2n=cos（pi*n/4）+cos（pi*n/8）;
X1k8=fft（x1n,8）;
X1k16=fft（x1n,16）;
X2k8=fft（x2n,8）;
X2k16=fft（x2n,16）;
subplot（2,2,1）;stem（w1,abs（X1k8）,'.'）;
title（'（1a）8点DFT[x1（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k8））]）;
subplot（2,2,2）;stem（w2,abs（X1k16）,'.'）;
title（'（1b）16点DFT[x1（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X1k16））]）;
subplot（2,2,3）;stem（w1,abs（X2k8）,'.'）;
title（'（2a）8点DFT[x2（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k8））]）;
subplot（2,2,4）;stem（w2,abs（X2k16）,'.'）;
title（'（2b）16点DFT[x2（n）]'）;
xlabel（'w/pi'）;ylabel（'幅度'）;
axis（[0,2,0,1.2*max（abs（X2k16））]）;
（3）对模拟周期信号进行谱分析
采样频率
，变换区间N=16，32，64三种情况进行谱分析，并画出幅频特性曲线。
clearall;closeall;
Fs=64;T=1/Fs;
N=16;n=0:
N-1;
xnT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;对x（t）取样；
Xk16=fft（xnT,16）;
Xk16=fftshift（Xk16）;将0频率移到频谱中心；
Tp=N*T;F=1/Tp;求频率分辨率F
k=-N/2:
N/2-1;fk=k*F;
subplot（3,1,1）;
stem（fk,abs（Xk16）,'.'）;
title（'（a）16点|DFT[x（nT）]|'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（Xk16））]）;
N=32;n=0:
N-1;
xnT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
Xk32=fft（xnT,32）;
Xk32=fftshift（Xk32）;
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot（3,1,2）;
stem（fk,abs（Xk32）,'.'）;
title（'（b）32点|DFT[x（nT）]|'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（Xk32））]）;
N=64;n=0:
N-1;
xnT=cos（8*pi*n*T）+cos（16*pi*n*T）+cos（20*pi*n*T）;
Xk64=fft（xnT,64）;
Xk64=fftshift（Xk64）;
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot（3,1,3）;
stem（fk,abs（Xk64）,'.'）;
title（'（c）64点|DFT[x（nT）]|'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;
axis（[-N*F/2-1,N*F/2-1,0,1.2*max（abs（Xk64））]）;
思考题：
1、对于周期序列，如果周期不知道，如何用FFT进行谱分析？
2、如何选择FFT的变换区间（包括非周期信号和周期信号）？
3、当N=8时，x2（n）和x3（n）的幅频特性会相同吗？
为什么？
N=16吗？
实验四：
IIR数字滤波器设计与软件实现
一、实验目的（略）
二、实验原理（略）
三、实验内容及步骤（略）
四、相关程序实现
1、调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号s（t），三路信号时域混叠无法分离，但频域是分离的，所以可以通过滤波实现信号的分离。
functionst=mstg编写程序mstg，供后面调用；
N=1600;Fs=10000;T=1/Fs;Tp=N*T;N为信号s（t）的长度；采样频率Fs=10kHz；
t=0:
T:
（N-1）*T;
k=0:
N-1;
f=k/Tp;
fc1=Fs/10;第一路调幅信号的载波频率fc1=1kHz;
fm1=fc1/10;第一路调幅信号的调制信号频率fm1=100Hz;
fc2=Fs/20;第二路调幅信号的载波频率fc2=500Hz;
fm2=fc2/10第二路调幅信号的调制信号频率fm2=50Hz;
fc3=Fs/40;第三路调幅信号的载波频率fc3=250Hz;
fm3=fc3/10;第三路调幅信号的调制信号频率fm3=25Hz;
xt1=cos（2*pi*fm1*t）.*cos（2*pi*fc1*t）;产生第一路调幅信号；
xt2=cos（2*pi*fm2*t）.*cos（2*pi*fc2*t）;产生第二路调幅信号；
xt3=cos（2*pi*fm3*t）.*cos（2*pi*fc3*t）;产生第三路调幅信号；
st=xt1+xt2+xt3;三路信号叠加形成时域混叠信号；
fxt=fft（st,N）;计算信号st的频谱
subplot（3,1,1）;
plot（t,st）;grid;plot线性x-y坐标绘图；图上加坐标网格；
xlabel（'t/s'）;ylabel（'s（t）'）;
axis（[0,Tp/2,min（st）,max（st）]）;
title（'（a）s（t）的波形'）;
subplot（3,1,3）;
stem（f,abs（fxt）/max（abs（fxt））,'.'）;grid;stem离散序列绘图；
title（'（b）s（t）的频谱'）;
axis（[0,Fs/5,0,1.2]）;
xlabel（'f/Hz'）;ylabel（'幅度'）;
由所产生的图形可以看出，频域内有三组谱线，载波频率分别是250Hz，500Hz和1000Hz，所以可以用低通滤波器滤出载波频率为250Hz的信号，用带通滤出载波频率为500Hz的信号，用高通滤出载波频率为1000Hz的信号。
2、
（1）要求将s（t）中三路调幅信号分离，通过观察s（t）的幅频特性曲线，分别确定可以分离s（t）中三路抑制载波单频调幅信号的三个滤波器（低通、带通、高通）的通带截止频率和阻带截止频率。
要求滤波器的通带最大衰减
为0.1dB，阻带最小衰减
为60dB。
（2）用MATLAB滤波器设计函数ellipord和ellip分别设计这三个椭圆滤波器，并绘制显示其幅频特性曲线。
（3）调用滤波器实现函数filter，用三个滤波器分别对s（t）进行滤波，分离出三路不同载波频率的调幅信号y1（n）、y2（n）、y3（n），并绘制三个信号的时域波形没观察分离效果。
由上个程序的运行结果可以看出载波频率分别是250Hz，500Hz和1000Hz，带宽分别为50Hz、100Hz、200Hz。
所以，三个滤波器的指标参数选取如下：
1对载波频率为250Hz的调幅信号，用低通分离，其指标为：
通带截止频率
，通带最大衰减
；
阻带截止频率
，阻带最小衰减
②对载波频率为500Hz的调幅信号，用带通分离，其指标为：
通带截止频率
，
，通带最大衰减
；
阻带截止频率
，
，阻带最小衰减
③对载波频率为1000Hz的调幅信号，用高通分离，其指标为：
通带截止频率
，通带最大衰减
；
阻带截止频率
，阻带最小衰减
说明：
为了使滤波器的阶数尽可能低，每个滤波器的边界频率的选择原则是尽量使滤波器的过渡带更宽一些。
clearall;closeall;
Fs=10000;T=1/Fs;M=1600;t=0:
T:
（M-1）*T;
st=mstg;调用子程序mstg；
==============以下为低通滤波器的设计与实现======================
fp=280;fs=450;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;获取数字滤波器的指标；
[N,wp]=ellipord（wp,ws,rp,rs）;计算椭圆DF的阶数N；
[B,A]=ellip（N,rp,rs,wp）;计算椭圆低通DF系统函数系数向量B和A；
[h1w,w]=freqz（B,A）;计算DF的频率响应
y1t=filter（B,A,st）;对信号st进行滤波
figure
（2）;subplot（3,1,1）;figure创建图形窗；
plot（w/pi,20*log10（abs（h1w）））;gridgrid图上加坐标网络；
axis（[0,1,-100,10]）;
xlabel（'w/pi'）;ylabel（'幅度（dB）'）;
subplot（3,1,3）;
plot（t,y1t）;grid
axis（[0,0.1,-1.5,1.5]）;
xlabel（'t/s'）;ylabel（'y1（t）'）;
===============以下为带通滤波器的设计与实现=========================
fpl=440;fpu=560;fsl=275;fsu=900;rp=0.1;rs=60;
wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];
[N,wp]=ellipord（wp,ws,rp,rs）;
[B,A]=ellip（N,rp,rs,wp）;
[h2w,w]=freqz（B,A）;
y2t=filter（B,A,st）;
figure（3）;subplot（3,1,1）;
plot（w/pi,20*log10（abs（h2w）））;grid
axis（[0,1,-100,10]）;
xlabel（'w/pi'）;ylabel（'幅度（dB）'）;
subplot（3,1,3）;
plot（t,y2t）;grid
axis（[0,0.1,-1.5,1.5]）;
xlabel（'t/s'）;ylabel（'y2（t）'）;
===============以下为高通滤波器的设计与实现====================
fp=890;fs=600;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
[N,wp]=ellipord（wp,ws,rp,rs）;
[B,A]=ellip（N,rp,rs,wp,'high'）;
[h3w,w]=freqz（B,A）;
y3t=filter（B,A,st）;
figure（4）;subplot（3,1,1）;
plot（w/pi,20*log10（abs（h3w）））;grid
axis（[0,1,-100,10]）;
xlabel（'w/pi'）;ylabel（'幅度（dB）'）;
subplot（3,1,3）;
plot（t,y3t）;grid
axis（[0,0.1,-1.5,1.5]）;
xlabel（'t/s'）;ylabel（'y3（t）'）;
升级会员 