北理工数字信号处理实验.docx
《北理工数字信号处理实验.docx》由会员分享，可在线阅读，更多相关《北理工数字信号处理实验.docx（45页珍藏版）》请在冰点文库上搜索。
本科实验报告
实验名称：
数字信号处理实验
课程名称：
数字信号处理实验
实验时间：
任课教师：
实验地点：
实验教师：
实验类型：
□原理验证
□综合设计
□自主创新
学生姓名：
学号/班级：
组号：
学院：
同组搭档：
专业：
成绩：
实验1利用DFT分析信号频谱
一、实验目的
1.加深对DFT原理的理解。
2.应用DFT分析信号频谱。
3.深刻理解利用DFT分析信号频谱的原理，分析现实过程现象及解决办法。
二、实验原理
1、DFT和DTFT的关系
有限长序列的离散时间傅里叶变换在频率区间的N个等分点上的N个取样值可以由下式表示：
由上式可知，序列的N点DFT，实际上就是序列的DTFT在N个等间隔频率点上样本。
2、利用DFT求DTFT
方法1：
由恢复出的方法如图2.1所示：
图2.1.由N点DFT恢复频谱DTFT的流程
由图2.1所示流程图可知：
由式2-2可以得到
其中为内插函数
方法2：
然而在实际MATLAB计算中，上诉插值公式不见得是最好的方法。
由于DFT是DTFT的取样值，其相邻的两个频率样本点的间距为,所以如果我们增加数据的长度N，使得得到的DFT谱线就更加精细，其包络就越接近DTFT的结果，这样可以利用DFT来近似计算DTFT。
如果没有更多的数据，可以通过补零来增加数据长度。
3、利用DFT分析连续时间信号的频谱
采用计算机分析连续时间信号的频谱，第一步就是把连续时间信号离散化，这里需要进行连个操作：
一是采样，二是截断。
对于连续非周期信号，按采样间隔T进行采样，截取长度为M，那么
对进行N点的频率采样，得到
因此，可以将利用DFT分析连续非周期信号频谱的步骤归纳如下：
（1）确定时域采样间隔T，得到离散序列；
（2）确定截取长度M，得到M点离散序列，这里的为窗函数。
（3）确定频域采样点数N，要求。
（4）利用FFT计算离散序列的N点DFT，得到。
（5）根据式（2-6）由计算采样点的近似值。
采用上诉方法计算的频谱，需要注意如下三点问题：
（1）频谱混叠。
如果不满足采样定理的条件，频谱会很出现混叠误差。
对于频谱无限宽的信号，应考虑覆盖大部分主要频率的范围。
（2）栅栏效应和频谱分辨率。
使用DFT计算频谱，得到的结果只是N个频谱样本值，样本值之间的频谱是未知的，就像通过一个栅栏观察频谱，称为“栅栏效应”。
频谱分辨率与记录长度成正比，提高频谱分辨率，就要增加记录时间。
（3）频谱泄露。
对于信号截断会把窗函数的频谱会引入到信号频谱中，造成频谱泄露。
解决这问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差。
因此，要合理选取采样间隔和截取长度，必要时还需考虑适当的窗。
对于连续周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上诉方法近似计算。
4、可能用到MATLAB函数与代码
实验中的DFT运算可以采用MATLAB中提供的FFT来实现。
DTFT可以利用MATLAB矩阵运算的方法进行计算。
三、实验题目
1.，完成如下要求：
（1）计算其DTFT，并画出区间的波形。
（2）计算4点DFT，并把结果显示在
（1）所画的图形中。
（3）对补零，计算64点DFT，并显示结果。
（4）是否可以由DFT计算DTFT，如果可以，请编程实现。
2.考察序列
（1）时，用DFT估计的频谱；将补零加长到长度为100点序列用DFT
估计的频谱。
要求画出相应波形。
（2）时，用DFT估计x（n）的频谱，并画出波形。
3.已知信号，其中，，。
从的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。
4.利用DFT近似分析连续时间信号xt=e-0.1tu（t）的频谱（幅度谱）。
分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定合适的参数。
四、实验代码、实验结果及分析
1.
（1）
>>n=0:
3;
>>x=[2-111];
>>w=-pi:
0.01*pi:
pi;
>>X=x*exp（-j*n'*w）;
>>subplot（211）;
>>plot（w,abs（X））;
>>xlabel（'\Omega/\pi'）;
>>title（'Magnitude'）;
>>axistight;
>>subplot（212）;
>>plot（w,angle（X））;
>>xlabel（'\Omega/\pi'）;
>>title（'Phase'）;
>>axistight;
（2）
>>Xk=fft（x）;
>>subplot（211）;
>>holdon;
>>stem（2*pi/4*n,abs（Xk）,'filled'）;
>>subplot（212）;
>>holdon;
>>stem（2*pi/4*n,angle（Xk）,'filled'）;
（3）
>>x=[2-111zeros（1,60）];
>>n=0:
63;
>>Xk=fft（x）;
>>subplot（211）;
>>stem（n,abs（Xk）,'filled'）;
>>xlabel（'\Omega/\pi'）;
>>title（'Magnitude'）;
>>axistight;
>>subplot（212）;
>>stem（n,angle（Xk）,'filled'）;
>>xlabel（'\Omega/\pi'）;
>>title（'Phase'）;
>>axistight;
（4）可以由DFT计算DTFT，只需采样点足够即可逼近DTFT，可采用增加补零个数的方法实现。
2.
（1）
>>n=0:
10;
>>x=cos（0.48*pi*n）+cos（0.52*pi*n）;
>>Xk=fft（x）;
>>subplot（211）;
>>stem（n,abs（Xk）,'filled'）;
>>xlabel（'\Omega'）;
>>title（'Magnitude'）;
>>n1=0:
99;
>>x1=[xzeros（1,89）];
>>Xk1=fft（x1）;
>>subplot（212）;
>>stem（n,abs（Xk1）,'filled'）;
>>xlabel（'\Omega'）;
>>title（'Magnitude'）;
>>n=0:
10;
>>x=cos（0.48*pi*n）+cos（0.52*pi*n）;
>>Xk1=fft（x）;
>>subplot（211）;
>>stem（n,angle（Xk）,'filled'）;
>>xlabel（'\Omega'）;
>>title（'Phase'）;
>>n1=0:
99;
>>x1=[xzeros（1,89）];
>>Xk1=fft（x1）;
>>subplot（212）;
>>stem（n1,angle（Xk1）,'filled'）;
>>xlabel（'\Omega'）;
>>title（'Phase'）;
（2）
（3）可通过增大采样点数，时域补零提高频谱分辨率。
3.
>>n=0:
20;
>>x=0.15*sin（2*pi*n）+sin（2*pi*2*n）-0.1*sin（2*pi*3*n）;
>>Xk=fft（x）;
>>subplot（211）;
>>stem（n,abs（Xk）,'filled'）;
>>xlabel（'\Omega'）;
>>title（'Magnitude'）;
>>subplot（212）;
>>stem（n,angle（Xk）,'filled'）;
>>xlabel（'\Omega'）;
>>title（'Phase'）;
4.
【修改采样间隔】
>>n=0:
10;
>>x=exp（-0.1*n）.*heaviside（n）;
>>Xk=fft（x）;
>>stem（n,abs（Xk）,'filled'）;
>>n=0:
50;
>>x=exp（-0.1*n）.*heaviside（n）;
>>Xk=fft（x）;
>>stem（n,abs（Xk）,'filled'）;
>>n=0:
100;
>>x=exp（-0.1*n）.*heaviside（n）;
>>Xk=fft（x）;
>>stem（n,abs（Xk）,'filled'）;
【修改截取长度】
>>n=0:
100;
>>x=exp（-0.1*n）.*（heaviside（n）-heaviside（n-50））;
>>Xk=fft（x）;
>>stem（n,abs（Xk）,'filled'）;
>>n=0:
100;
>>x=exp（-0.1*n）.*（heaviside（n）-heaviside（n-10））;
>>Xk=fft（x）;
>>stem（n,abs（Xk）,'filled'）;
五、实验总结
第一次数字信号处理实验，我了解了如何在MATLAB中使用DFT分析信号的频谱，对DFT采样、截取的概念和操作结果有了更直观的认识。
此外，使用MATLAB的FFT大大简化了DFT的计算，让我认识到了FFT算法的快捷。
实验2利用FFT计算线性卷积
一、实验目的
1.掌握利用FFT计算线性卷积的原理及具体实现方法。
2.加深理解重叠相加法和重叠保留法。
3.考察利用FFT计算线性卷积各种方法的适用范围。
二、实验原理
1.线性卷积与圆周卷积
设x（n）为L点序列，h（n）为M点序列，x（n）和h（n）的线性卷积为
（3-1）
的长度为L+M-1
x（n）和h（n）的圆周卷积为
（3-2）
圆周卷积与线性卷积相等而不产生交叠的必要条件为
N≥L+M+1（3-3）
圆周卷积定理：
根据DFT性质，x（n）和h（n）的N点圆周卷积的DFT等于它们的DFT的乘积：
（3-4）
2.快速卷积
快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。
可将快速卷积运算的步骤归纳如下：
（1）必须选择；为了能使用基-2算法，要求。
采用补零的办法使得x（n）和h（n）的长度均为N。
（2）计算x（n）和h（n）的N点FFT。
（3）组成乘积
（4）利用IFFT计算Y（k）的IDFT，得到线性卷积y（n）
3.分段卷积
我们考察单位取样响应为h（n）的线性系统，输入为x（n），输出为y（n），则
yn=xn*h（n）
当输入序列x（n）极长时，如果要等x（n）全部集齐时再开始进行卷积，会使输出有较大延时；如果序列太长，需要大量存储单元。
为此，我们把x（n）分段，为别求出每段的卷积，合在一起得到最后的总输出。
这称为分段卷积。
分段卷积可以细分为重叠保留法和重叠相加法。
重叠保留法：
设x（n）的长度为，h（n）的长度为M。
把序列x（n）分成多段N点序列，每段雨前一段重写M-1个样本。
并在第一个输入段前面补M-1个零。
计算每一段与h（n）的圆周卷积，其结果中前M-1个不等与线性卷积，应当舍去，只保留后面N-M+1个正确的输出样本，把它们合起来得到总的输出。
利用FFT实现重叠保留法的步骤如下：
（1）在x（n）前面填充M-1个零，扩大以后的序列为
（2）将x（n）分为若干段N点子段，设L=N-M+1为每一段的有效长度，则第i段的数据为：
（3）计算每一段与h（n）的N点圆周卷积，利用FFT计算圆周卷积
（4）舍去每一段卷积结果的前M-1个样本，连接剩下的样本得到卷积结果y（n）。
重叠相加法：
设h（n）长度为M，将信号x（n）分解成长为L的子段。
以表示没断信号，则：
每一段卷积的长度为L+M-1，所以在做求和时，相邻两段序列由M-1个样本重叠，即前一段的最后M-1个样本和下一段前M-1个样本序列重叠，这个重叠部分相加，再与不重叠的部分共同组成y（n）。
利用FFT实现重叠保留法的步骤如下：
（1）将x（n）分为若干L点子段。
（2）计算每一段与h（n）的卷积，根据快速卷积法利用FFT计算卷积。
（3）将各段相加，得到输出y（n）。
4、可能得到的MATLAB函数
实验中FFT运算可采用MATLAB中提供的函数fft来实现。
三、实验题目
假设要计算序列x（n）=u（n）-u（n-L）,0≤n≤L和h（n）=cos（0.2πn）,0≤n≤M的线性卷积完成以下实验内容。
1.设L=M,根据线性卷积的表达式和快速卷积的原理分别编程实现计算两个序列线性卷积的方法,比较当序列长度分别为8,16,32,64,256,512,1024时两种方法计算线性卷积所需时间。
2．当L=2048且M=256时比较直接计算线性卷积和快速卷积所需的时间,进一步考察当L=4096且M=256时两种算法所需的时间。
3.编程实现利用重叠相加法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与2题的结果进行比较。
4.编程实现利用重叠保留法计算两个序列的线性卷积,考察L=2048且M=256时计算线性卷积的时间,与2题的结果进行比较。
四、实验代码、实验结果及分析
1.
【线性卷积】
>>fori=1:
7
L=input（'L:
'）;
n=0:
L;
x=heaviside（n）-heaviside（n-L）;
h=cos（0.2*pi*n）;
tic
y=conv（x,h）;
toc
end
L:
8
时间已过0.077287秒。
L:
16
时间已过0.032760秒。
L:
32
时间已过0.000095秒。
L:
64
时间已过0.000060秒。
L:
256
时间已过0.000131秒。
L:
512
时间已过0.000230秒。
L:
1024
时间已过0.000846秒。
【快速卷积】
>>fori=1:
7
L=input（'L:
'）;
n=0:
L;
x=heaviside（n）-heaviside（n-L）;
h=cos（0.2*pi*n）;
tic
Xk=fft（x,L+1）;
Hk=fft（h,L+1）;
Yk=Xk.*Hk;
y=ifft（Yk）;
toc
end
L:
8
时间已过0.002417秒。
L:
16
时间已过0.001108秒。
L:
32
时间已过0.000384秒。
L:
64
时间已过0.000616秒。
L:
256
时间已过0.002367秒。
L:
512
时间已过0.000702秒。
L:
1024
时间已过0.000967秒。
在序列长度较短时，快速卷积优势明显，序列加长后出现了比线性卷积慢的情况
2.
线性卷积：
>>fori=1:
2
L=input（'L:
'）;
n1=0:
L;
M=input（'M:
'）;
n2=0:
M;
x=heaviside（n1）-heaviside（n1-L）;
h=cos（0.2*pi*n2）;
tic
y=conv（x,h）;
toc
end
L:
2048
M:
256
时间已过0.000533秒。
L:
4096
M:
256
时间已过0.000712秒。
快速卷积：
fori=1:
2
L=input（'L:
'）;
n1=0:
L;
M=input（'M:
'）;
n2=0:
M;
x=heaviside（n1）-heaviside（n1-L）;
h=cos（0.2*pi*n2）;
tic
Xk=fft（x,M）;
Hk=fft（h,M）;
Yk=Xk.*Hk;
y=ifft（Yk）;
toc
end
L:
2048
M:
256
时间已过0.000592秒。
L:
4096
M:
256
时间已过0.001827秒。
可见，当序列加长，快速卷积速度开始优于线性卷积
3.
>>tic
>>N=512;
>>m=0:
256;
>>h=cos（0.2*pi*m）;
>>n=0:
2048;
>>x=heaviside（n）-heaviside（n-2048）;
>>Lx=length（x）;
>>Lm=length（h）;
>>M=Lm-1;
>>L=N-M;
>>h=fft（h,N）;
>>K=ceil（Lx/L）;
>>fori=Lx:
K*L-1
x（i+1）=0;
end
>>Y=zeros（K,N）;
>>Y2=zeros（1,（K-1）*L+N）;
>>fork=0:
K-1
Xk=[x（k*L+1:
k*L+L）,zeros（1,M）];
Y（k+1,:
）=（ifft（fft（Xk）.*h））;
Y2（k*L+1:
k*L+N）=Y2（k*L+1:
k*L+N）+Y（k+1,:
）;
end
>>toc
时间已过0.000575秒。
可见，重叠相加法的速度相比第二题更快
4.
>>tic
>>N=512;
>>m=0:
256;
>>h=cos（0.2*pi*m）;
>>n=0:
2048;
>>x=heaviside（n）-heaviside（n-2048）;
>>Lx=length（x）;
>>Lm=length（h）;
>>M=Lm-1;
>>L=N-M;
>>h=fft（h,N）;
>>K=floor（（Lx+M-1）/L）+1;
>>z=（K）*L-Lx;
>>x1=[zeros（1,M）,x,zeros（1,z）];
>>Y=zeros（K,N）;
>>fork=0:
K-1
Xk=fft（x1（k*L+1:
k*L+N））;
Y（k+1,:
）=（ifft（Xk）.*h）;
end
>>Z=reshape（Y（:
M:
N）',1,[]）;
>>toc
时间已过0.000452秒。
可见，重叠保留法的速度相比第二题也更快
五、实验总结
这次FFT的实验让我对FFT的算法效率有了新的认识，并不是任何情况下FFT卷积的速度都优于直接计算线性卷积，当序列长度很长时反而会比线性卷积慢。
此外，我还对两种新的FFT方法：
重叠相加法、重叠保留法及其MATLAB实现有了基本的了解。
实验3IIR数字滤波器设计
一、实验目的
1.掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。
2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。
3.掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。
二、实验题目
1.设采样频率为fs=4kHz，采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率为fc=1kHz。
2.设采样频率为fs=10kHz，设计数字低通滤波器，满足如下指标
通带截止频率：
fp=1kHz，通带波动：
Rp=1dB
阻带截止频率：
fst=1.5kHz，阻带衰减：
As=15dB
要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。
结合实验结果，分别讨论采用上述设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。
三、实验代码、实验结果及分析
1.巴特沃斯模拟原型滤波器+脉冲响应不变法/双线性变换法
>>fp=1000;
>>fs=1500;
>>Rp=1;
>>As=15;
>>f=10000;
>>Wp=2*pi*fp;
>>Ws=2*pi*fs;
>>wp=2*f*tan（2*pi*fp/（2*f））;
>>ws=2*f*tan（2*pi*fs/（2*f））;
>>N1=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（Wp/Ws）））;
>>N2=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（wp/ws）））;
>>OmegaC1=Wp/（（10^（Rp/10）-1）^（1/（2*N1）））;
>>OmegaC2=wp/（（10^（Rp/10）-1）^（1/（2*N1）））;
>>[z,p,k]=buttap（N1）;
>>b1=k*OmegaC1^N1;
>>a1=poly（p*OmegaC1）;
>>[H,w]=freqs（b1,a1）;
>>subplot（221）;
>>plot（w/pi,abs（H））;
>>gridon;
>>axistight;
>>xlabel（'\Omega（\pi）'）;
>>ylabel（'|H（\Omega）|'）;
>>subplot（223）;
>>plot（w/pi,angle（H））;
>>gridon;
>>axistight;
>>xlabel（'\Omega（\pi）'）;
>>ylabel（'Phaseof|H（\Omega）|（\pi）'）;
>>[b2,a2]=butter（N2,OmegaC2,'S'）;
>>[H,w]=freqs（b2,a2）;
>>subplot（222）;
>>plot（w/pi,abs（H））;
>>gridon;
>>axistight;
>>xlabel（'\Omega（\pi）'）;
>>ylabel（'|H（\Omega）|'）;
>>subplot（224）;
>>plot（w/pi,angle（H））;
>>gridon;
>>axistight;
>>xlabel（'\Omega（\pi）'）;
>>ylabel（'Phaseof|H（\Omega）|（\pi）'）;
2.切比雪夫I型+脉冲响应不变法/双线性变换法
>>wp=0.2*pi;
>>ws=0.3*pi;
>>Rp=1;
>>As=15;
>>f=10000;
>>T=1/f;
>>e=（10^（Rp/10）-1）^（1/2）;
>>A=10^（As/20）;
>>wp1=wp/T;
>>ws1=ws/T;
>>N=ceil（（acosh（（A^2-1）^（1/2）/e））/（acos
升级会员 