数字信号处理实验报告.docx
《数字信号处理实验报告.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验报告.docx（33页珍藏版）》请在冰豆网上搜索。
数字信号处理实验报告
《数字信号处理》实验报告
课程名称
数字信号处理
学生姓名
指导教师
李宏
学院
信息科学与工程学院
专业班级
学号
2017年5月
实验一离散时间信号和系统响应
一.实验目的
1.掌握连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解
2.掌握求系统响应的方法
3.掌握时域离散系统的时域特性
4.利用卷积方法观察分析系统的时域特性
二.实验原理与方法
采样是连续信号数字化处理的第一个环节。
对采样过程的研究不仅可以了解采样前后信号时域和频域特性的变化以及信号信息不丢失的条件，而且可以加深对离散傅里叶变换、Z变换和序列傅里叶变换之间关系式的理解。
对连续信号
以T为采样间隔进行时域等间隔理想采样，形成采样信号
：
式中
为周期冲激脉冲，
的傅里叶变换
为：
上式表明将连续信号xa（t）采样后其频谱将变为周期的，周期为Ωs=2π/T。
也即采样信号的频谱是原连续信号xa（t）的频谱Xa（jΩ）在频率轴上以Ωs为周期，周期延拓而成的。
在时域中，描述系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。
一个时域离散线性时不变系统的输出与输入间的关系为：
这里，
为系统的输出序列，
为输入序列。
可以是无限长，也可以是有限长。
三.实验内容
1.时域采样定理的验证
给定模拟信号：
，
式中
。
对其进行采样，可得到采样序列
。
选择三种采样频率Fs=1kHz,300Hz,200Hz,观测时间选Tp=64ms，分别得到三个采样序列
。
编写程序计算这三个序列的幅度特性，并绘图显示其幅频特性曲线。
观察分析频谱混叠现象。
2.给定一个低通滤波器的差分方程为：
输入信号
a.分别求出
和
的系统响应，并画出其波形
b.求出系统的单位脉冲响应，画出其波形
3.给定系统的单位脉冲响应为
用线性卷积法求
分别对系统
和
的输出响应，并画出波形
A=444.128;
a=50*sqrt
（2）*pi;
w=50*sqrt
（2）*pi;
n=0:
63;
fs=1000;
c=A*exp（（-a）*n/fs）.*sin（w*n/fs）;
subplot（3,2,1）;
stem（n,c,'.'）;
xlabel（'n'）;
ylabel（'xa（n）'）;
title（'xa（n）的时域序列'）;
N=64;
k=0:
length（c）-1;
w=k*pi/100;
X=fft（c,N）;
subplot（3,2,2）;
plot（w/pi,abs（X））;
xlabel（'w/pi'）;
ylabel（'|X（jw）|'）;
title（'xa（n）的傅氏变换|X（jw）|'）;
A=444.128;
a=50*sqrt
（2）*pi;
w=50*sqrt
（2）*pi;
n=0:
63;
fs=300;
c=A*exp（（-a）*n/fs）.*sin（w*n/fs）;
subplot（3,2,3）;
stem（n,c,'.'）;
xlabel（'n'）;
ylabel（'xa（n）'）;
title（'xa（n）的时域序列'）;
N=64;
k=0:
length（c）-1;
w=k*pi/100;
X=fft（c,N）;
subplot（3,2,4）;
plot（w/pi,abs（X））;
xlabel（'w/pi'）;
ylabel（'|X（jw）|'）;
title（'xa（n）的傅氏变换|X（jw）|'）;
A=444.128;
a=50*sqrt
（2）*pi;
w=50*sqrt
（2）*pi;
n=0:
63;
fs=200;
c=A*exp（（-a）*n/fs）.*sin（w*n/fs）;
subplot（3,2,5）;
stem（n,c,'.'）;
xlabel（'n'）;
ylabel（'xa（n）'）;
title（'xa（n）的时域序列'）;
N=64;
k=0:
length（c）-1;
w=k*pi/100;
X=fft（c,N）;
subplot（3,2,6）;
plot（w/pi,abs（X））;
xlabel（'w/pi'）;
ylabel（'|X（jw）|'）;
title（'xa（n）的傅氏变换|X（jw）|'）;
内容一：
调用filter解差分方程，由系统对u（n）的响应判断稳定性
A=[1,-0.9];
B=[0.05,0.05];
x1n=[11111111zeros（1,122）];
x2n=ones（1,122）;
hn=impz（B,A,58）;
subplot（2,2,1）;
y='hn';
stem（hn,'r','.'）;
xlabel（'n'）;
ylabel（'h（n）'）;
title（'（a）系统单位脉冲响应h（n）'）;
y1n=filter（B,A,x1n）;
subplot（2,2,2）;
y='y1n';
stem（y1n,'g','.'）;
xlabel（'n'）;
ylabel（'h（n）'）;
title（'（b）系统对R8（n）的响应y1（n）'）;
y2n=filter（B,A,x2n）;
subplot（2,2,4）;
y='y1n';
stem（y2n,'b','.'）;
xlabel（'n'）;
ylabel（'y2（n）'）;
title（'（c）系统对u（n）的响应y2（n）'）;
内容一:
系统的单位冲响应的波形如下图（a）所示，系统和的响应序列的波形如下图（b）和图（c）
%内容2:
调用conv函数计算卷积
x1n=[11111111];
h1n=[ones（1,10）zeros（1,22）];
h2n=[12.52.51zeros（1,22）];
y21n=conv（h1n,x1n）;
y22n=conv（h2n,x1n）;
figure
（2）;
subplot（2,2,1）;
y='h1（n）';
stem（h1n,'b','.'）;
axis（[02201.2]）;
title（'（d）系统单位脉冲响应h1（n）'）;
subplot（2,2,2）;
y='y21n';
stem（y21n,'r','.'）;
axis（[05208.2]）;
title（'（e）h1（n）与R8（n）的卷积y21（n）'）;
subplot（2,2,3）;
y='h2（n）';
stem（h2n,'b','.'）;
axis（[02203.2]）;
title（'（f）系统对单位脉冲响应h2（n）'）;
subplot（2,2,4）;
y='y22n';
stem（y22n,'r','.'）;
gridon;
title（'（g）h2（n）与R8（n）的卷积y22（n）'）;
四.实验思考
1.在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？
它们所对应的模拟频率是否相同？
为什么？
；
对于生活中的连续信号（图像，语音等），采样是十分重要的，采样频率必须满足采样定理（f>2fs）才能安全重构原来的信号，当采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是不同的，它们所对应的模拟频率也是不同的。
因为采样的频率不同，其傅里叶变换所能包含的范围也不同。
而我们是在抽样的过程中考虑采样频率，在FFT中其实可以不考虑，归一化都可以，所以说提高采样率我们的重点是能不能够重构原来的信号，如果原连续信号的最高频率分量或者带宽很宽的话，那么按照采样定理就要提高采样率了。
2.如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应？
如何求？
直接进行线性卷积是不可行的，因为输入信号为无限长序列。
但在实验中函数fftt的调用格式为y=fftfilt（b,x）该格式是利用基于FFT的重叠相加法对数据进行滤波,这种频域滤波技术只对FIR滤波器有效.该函数是通过向量b描述的滤波器对x数据进行滤波.
x是等待滤波的信号；
b是FIR滤波器的H（z）的分子多项式系数
3.如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化？
用前面第二个实验结果进行分析说明
如果信号经过低通滤波器,把信号的高频分量滤掉,时域信号将变得更加平滑。
并且由实验内容
（2）可见,经过系统低通滤波使输入信号、和的阶跃变化变得缓慢上升与下降.
实验总结
（1）时域求系统响应的方法有两种，第一种是通过解差分方程求得系统的响应输出，注意合理的选择系统的初始条件，第二种是已知系统的单位冲激响应，通过求系统输入与单位脉冲响应的线性卷积求得系统的输出。
（2）实验中要检查系统的稳定性，其方法是在输入端加上单位阶跃序列，观察系统的输出波形，如果波形的幅值稳定在一个常数值附近，则系统稳定，否则系统不稳定。
上面第三个实验是稳定的。
实验二用FFT对信号作频谱分析
一.实验目的
1.进一步加深DFT算法原理和基本性质的理解
2.掌握用FFT对连续信号和时域离散信号进行频谱分析的方法
3.了解用FFT进行频谱分析时可能出现的分析误差及其原因，以便正确应用FFT。
二.实验原理
（补充离散傅里叶变换的性质）
用FFT对信号作频谱分析是学习数字信号处理的重要内容，经常需要进行谱分析的信号是模拟信号和时域离散信号。
对信号进行谱分析的重要问题是频谱分辨率F和分析误差。
频谱分辨率直接和FFT的变换区间N有关，FFT能够实现的频率分辨率是2/N，因此要求2/NF。
可以根据此式选择FFT的变换区间N。
误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时，离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。
周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。
如果不知道信号周期，可以尽量选择信号的观察世间长一些。
对模拟信号进行谱分析时，首先要按照采样定理将其变为时域离散信号。
如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。
三.实验内容
1.对以下给出的各序列进行谱分析：
选择FFT的变换区间N为8和16两种情况进行频谱分析。
2.对以下各周期序列进行频谱分析
选择FFT的变换区间N为8和16两种情况分别对以上序列进行频谱分析。
3.对模拟周期信号进行频谱分析
选择采样频率Fs=64Hz，对变换区间N=16，32，64三种情况进行谱分析
1.对以下序列进行谱分析:
选择FFT的变换区间N为8和16两种情况进行频谱分析。
分别打印其幅频特性曲线,并进行对比、分析和讨论。
代码如下：
x1n=[ones（1,4）];%产生R4（n）序列向量
X1k8=fft（x1n,8）;
X1k16=fft（x1n,16）;
N=8;
f=2/N*（0:
N-1）;
figure
（1）;
subplot（1,2,1）;
stem（f,abs（X1k8）,'.'）;
title（'（1a）8点DFT[x_1（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
N=16;
f=2/N*（0:
N-1）;
figure
（1）;
subplot（1,2,2）;
stem（f,abs（X1k16）,'.'）;
title（'（1a）16点DFT[x_1（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
M=8;
xa=1:
（M/2）;
xb=（M/2）:
-1:
1;
x2n=[xa,xb];
x3n=[xb,xa];
X2k8=fft（x2n,8）;
X2k16=fft（x2n,16）;
X3k8=fft（x3n,8）;
X3k16=fft（x3n,16）;
figure
（2）;
N=8;
f=2/N*（0:
N-1）;
subplot（2,2,1）;
stem（f,abs（X2k8）,'.'）;
title（'（2a）8点DFT[x_1（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
subplot（2,2,3）;
stem（f,abs（X3k8）,'.'）;
title（'（3a）8点DFT[x_1（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
N=16;
f=2/N*（0:
N-1）;
subplot（2,2,2）;
stem（f,abs（X2k16）,'.'）;
title（'（2a）16点DFT[x_1（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
subplot（2,2,4）;
stem（f,abs（X3k16）,'.'）;
title（'（3a）16点DFT[x_1（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
N=8;
n=0:
N-1;
x4n=cos（pi*n/4）;
x5n=cos（pi*n/4）+cos（pi*n/8）;
X4k8=fft（x4n,8）;
X4k16=fft（x4n,16）;
X5k8=fft（x5n,8）;
X5k16=fft（x5n,16）;
figure（3）;
N=8;
f=2/N*（0:
N-1）;
subplot（2,2,1）;
stem（f,abs（X4k8）,'.'）;
title（'（4a）8点DFT[x_4（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
subplot（2,2,3）;
stem（f,abs（X5k8）,'.'）;
title（'（5a）16点DFT[x_5（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
N=16;
f=2/N*（0:
N-1）;
subplot（2,2,2）;
stem（f,abs（X4k16）,'.'）;
title（'（4a）16点DFT[x_4（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
subplot（2,2,4）;
stem（f,abs（X5k16）,'.'）;
title（'（5a）16点DFT[x_5（n）]'）;
xlabel（'w/pai'）;
ylabel（'幅度'）;
四.实验思考
1.在N=8时，
和
的幅频特性会相同吗？
为什么？
N=16时呢？
答：
当N=8时，幅频特性相同，由时域循环移位定理，x2（n）和x3（n）满足时域循环移位定理。
当N=16时，它们的幅频特性不会相同，因为，x2（n）和x3（n）后都要补零，x2（n）进行循环移位后的序列跟x3（n）不相同。
2.对于周期序列，如果周期不知道，如何用FFT进行谱分析？
答：
周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差要求就继续将截取长度加倍,重复比较,直到结果满足要求
3.如何选择FFT的变换区间？
（包括非周期信号和周期信号）
答：
一、对于非周期信号：
有频谱分辨率F，而频谱分辨率直接和FFT的变换区间有关，因为FFT能够实现的频率分辨率是2π/N...因此有最小的N>2π/F。
就可以根据此式选择FFT的变换区间。
二、对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。
答：
周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。
如果不知道信号周期，可以尽量选择信号的观察时间长一些。
对模拟信号进行频谱分析时，首先要按照采样定理将其变成时域离散信号。
如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的普分析进行。
实验三IIR数字滤波器的设计
一.实验目的
1.熟悉用双线性变换法设计IIR数字滤波器的原理和方法
2.掌握数字滤波器的Matlab实现方法
3.通过观察对实际心电图信号的滤波作用，获得数字滤波的感性认识
二.实验原理
设计IIR数字滤波器一般采用间接设计法——脉冲响应不变法和双线性变换法，应用最广泛的是双线性变换法。
脉冲响应不变法的基本思想是：
使数字滤波器的单位脉冲响应h（n）近似于模拟滤波器的单位脉冲响应ha（t），即使
其S平面和Z平面的映射关系为：
双线性变换法的基本思想是：
使描述数字滤波器的差分方程近似描述模拟滤波器的微分方程
S平面和Z平面的映射关系为：
双线性变换法中的频率变换是一种非线性变换，这种非线性引起的幅频特性畸变可通过预变形矫正法而得到校正。
设计IIR数字滤波器的一般步骤:
（1）确定所需类型数字滤波器的技术指标：
通带截止频率ωp、通带衰减αp、阻带截止频率ωs、阻带衰减αs。
（2）将所需类型数字滤波器的技术指标转换成相应类型模拟滤波器的技术指标。
（3）设计该类型模拟滤波器
（4）通过复频率变换将模拟滤波器转换成所需类型的数字滤波器。
三.实验内容
1.分别用脉冲响应不变法和双线性变换法设计一个巴特沃斯低通IIR数字滤波器，设计指标参数为：
在通带内频率低于0.2时，最大衰减小于1dB，在阻带内[0.3，]频率区间上，最小衰减大于15dB。
观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。
比较这两种方法的优缺点。
Fp=0.1;
Fs=0.15;
Ft=1;
Nn=256;
As=15;
Ap=1
wp=2*pi*Fp;
ws=2*pi*Fs;
[n,wn]=buttord（wp,ws,Ap,As,'s'）;%求低通滤波器的阶数和截止频率；n为滤波器最小阶数，wn为其截止频率，
%wp，ws分别为通带频率和截止频率，Ap为通带最大衰减，As阻带最小衰减。
[b,a]=butter（n,wn,'s'）;%求S域的频率响应的参数
[bz,az]=impinvar（b,a,Ft）;%脉冲响应不变法实现S域到Z域的变换
[h,f]=freqz（bz,az,Nn,Fs）;%根据参数求出频率响应
figure
（1）
subplot（2,1,1）;plot（f,20*log10（abs（h）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;
gridon;
subplot（2,1,2）;plot（f,180/pi*unwrap（angle（h）））;
xlabel（'频率/Hz'）;ylabel（'相位/^o'）;
gridon;
title（'脉冲响应不变法'）;
[bz,az]=bilinear（b,a,Ft）;%双线性变换实现S域到Z域的变换
[h,f]=freqz（bz,az,Nn,Fs）;%根据参数求出频率响应
figure
（2）
subplot（2,1,1）;plot（f,20*log10（abs（h）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;
gridon;
subplot（2,1,2）;plot（f,180/pi*unwrap（angle（h）））;
xlabel（'频率/Hz'）;ylabel（'相位/^o'）;
gridon;
title（'双线性变换法'）;
3.心电图
xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-19,-38,-60,84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
n=0:
55;
Fp=0.1;
Fs=0.15;
Ft=1;
Nn=256;
As=15;
Ap=1
wp=2*pi*Fp;
ws=2*pi*Fs;
[k,wk]=buttord（wp,ws,Ap,As,'s'）;%求低通滤波器的阶数和截止频率；k滤波器最小阶数，wk截止频率，
%wp，ws分别为通带频率和截止频率，Ap为通带最大衰减，As阻带最小衰减。
[b,a]=butter（k,wk,'s'）;%求S域的频率响应的参数
[bz,az]=bilinear（b,a,Ft）;%双线性变换实现S域到Z域的变换
[h,f]=freqz（bz,az,Nn,Fs）;%根据参数求出频率响应
figure
（1）
stem（n,xn）;
title（'原始心电信号'）;
y1=filter（bz,az,xn）;
figure
（2）
plot（n,y1）;
title（'滤波后的心电信号'）
四.实验思考
1.用双线性变换法设计数字滤波器过程中，变换公式
中T的取值，对设计结果有无影响？
为什么？
答：
无影响，因为数字滤波器的传输函数H（ejw）以2
为周期，最高频率在
处，因此，
<
，按照线性关系
那么一定满足
因此T可以任选。
在数字滤波器的设计过程中，设计问题起始于数字滤波器上的传输函数
的技术要求，当这些技术要求通过变换公式
得到模拟滤波器的传输函数HC（s）的技术要求并设计出模拟滤波，然后模拟滤波器的传输函数HC（s）再通过公式变换得到数字滤波器的传输函数
，这样取值T对设计的影响就被抵消，因此T的取值对设计结果没有影响。
实验四FIR数字滤波器的设计
一.实验目的
1.掌握用窗函数法设计FIR数字滤波器的原理和方法
2.熟悉线性相位FIR数字滤波器特性
3.了解各种窗函数对滤波特性的影响
二.实验原理
窗函数法设计FIR滤波器的步骤为：
（1）构建希望逼近的理想频率响应函数
及技术指标
（2）求滤波器的单位脉冲响应
如果
复杂，可对
从
采样M个点，采样值为
，则：
（3）根据对过渡带及阻带衰减的要求，选择窗函数的形式，并估计窗口宽度N，设要求的过渡带宽为
，则
（4）计算滤波器的单位脉冲响应
（5）求H（ejω），分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。
窗函数傅里叶变换W（ejω）的主瓣决定了H（ejω）过渡带宽，W（ejω）的旁瓣大小和多少决定了H（ejω）在通带和阻带范围内波动幅度，常用的几种窗函数有：
矩形窗；Hanning窗；Hamming窗；Blackmen窗；Kaiser窗
三.实验内容
1.用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率
。
窗口长度N=15，33。
要求在两种窗口长度情况下，分别求出h（n），打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽，总结窗
升级会员 