北理工数字信号处理实验报告.docx

北理工数字信号处理实验报告.docx

《北理工数字信号处理实验报告.docx》由会员分享，可在线阅读，更多相关《北理工数字信号处理实验报告.docx（37页珍藏版）》请在冰豆网上搜索。

本科实验报告

实验名称：

数字信号处理实验

课程名称：

数字信号处理实验

实验时间：

周二8~10节

任课教师：

范哲意

实验地点：

4—423

实验教师：

范哲意

实验类型：

□原理验证

□综合设计

□自主创新

学生姓名：

任科飞

学号/班级：

1120121159/05111251

组号：

学院：

信息与电子

同组搭档：

专业：

信息工程

成绩：

实验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）的频谱（幅度谱）。

分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定合适的参数。

四、实验代码及实验结果

第一题

实验代码

x=[2,-1,1,1];

n=0:

3;

n1=0:

63;

z=zeros（1,60）;

x1=[xz];

y=fft（x）;

y1=fft（x1）;

w=-pi:

0.01*pi:

pi;

X=x*exp（-j*n'*w）;

subplot（211）;

plot（w,abs（X））;

xlabel（'\Omega/\pi'）;

title（'Magnitude'）;

holdon;

stem（n,y,'filled'）;

axistight;

subplot（212）;

plot（w,angle（X）/pi）;

xlabel（'\Omega/\pi'）;

title（'Phase'）;

axistight;

figure;

stem（n1,y1,'filled'）;

实验结果

分析

可以由DFT计算DTFT。

增加补零的个数，提高采样密度，可以由DFT近似计算DTFT。

第二题

实验代码

n=0:

10;

x=cos（0.48*pi*n）+cos（0.52*pi*n）;

n1=0:

99;

n2=0:

100;

z=zeros（1,89）;

x1=[xz];

x2=cos（0.48*pi*n2）+cos（0.52*pi*n2）;

y=fft（x）;

y1=fft（x1）;

y2=fft（x2）;

subplot（211）;

stem（n,y,'filled'）;

title（'0<=n<=10'）;

subplot（212）;

stem（n1,y1,'filled'）;

title（'0<=n<99'）;

figure;

stem（n2,y2,'filled'）;

title（'0<=n<=100'）;

实验结果

（3）分析

可以通过增大截取长度和增加补零的个数来提高频谱分辨率，但是补零不能够增加分辨力。

第三题

实验代码

>>n=0:

39;

>>x1=0.15*sin（2*pi*n）+sin（4*pi*n）-0.1*sin（6*pi*n）;

>>X1=fft（x1,40）;

>>subplot（311）

>>stem（n,abs（X1））;

>>x2=0.15*sin（2*pi*（1/6）*n）+sin（4*pi*（1/6）*n）-0.1*sin（6*pi*（1/6）*n）;

>>X2=fft（x2,40）;

>>subplot（312）

>>stem（n,abs（X2））;

>>x3=0.15*sin（2*pi*0.1*n）+sin（4*pi*0.1*n）-0.1*sin（6*pi*0.1*n）;

>>X3=fft（x3,40）;

>>subplot（313）

>>stem（n,abs（X3））;

实验结果

分析：

减小采样周期，可以使频谱分辨率提高。

第四题

实验代码

n=0:

100;

x=exp（-0.1*n）;

X1=fft（x）;

k=0:

100;

stem（k,abs（X1））;

xlabel（'k'）;

实验结果

实验代码

n=0:

2:

100;

x=exp（-0.1*n）;

X2=fft（x）;

k=0:

50;

stem（k,abs（X2））;

xlabel（'k'）;

实验结果

实验代码

>>n=0:

50;

x=exp（-0.1*n）;

X3=fft（x）;

k=0:

50;

stem（k,abs（X3））;

xlabel（'k'）;

实验结果

实验代码

>>n=0:

2:

50;

x=exp（-0.1*n）;

X4=fft（x）;

k=0:

25;

stem（k,abs（X4））;

xlabel（'k'）;

实验结果

实验代码

>>n=0:

150;

x=exp（-0.1*n）;

X5=fft（x）;

k=0:

150;

stem（k,abs（X5））;

xlabel（'k'）;

实验结果

实验代码

>>n=0:

2:

150;

x=exp（-0.1*n）;

X5=fft（x）;

k=0:

75;

stem（k,abs（X5））;xlabel（'k'）;

实验结果

分析：

采样区间可以为0到100，采样间隔为2。

五、心得与体会

通过本次实验，我们掌握并加深对DFT原理的理解并且学会应用DFT分析信号频谱，在此基础上利用DFT分析信号频谱的原理，掌握了分析现实现象的步骤及办法。

并且，通过这次的实验，我对DFT的原理确实有了更深的了解，对DFT和DTFT之间的关系也有了更多的认识：

序列x（n）的N点DFTX（K），实际上就是x（n）序列的DTFT在N个等间隔频率点上的样本X（K）。

所以，我们可以通过增加数据的长度N或通过补零来使DFT更加接近DTFT的结果。

这样就可以利用DFT近似计算DTFT。

实验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题的结果进行比较。

四、实验代码及实验结果

第一题

实验代码

L=input（'length='）;

n1=0:

2*L-1;

n2=0:

L-1;

x=ones（1,L）;

h=cos（0.2*pi*n2）;

z=zeros（1,L-1）;

x1=[xz];

h1=[hz];

tic;

y1=conv（x,h）;

toc;

tic;

X=fft（x1）;

H=fft（h1）;

Y=X.*H;

y=ifft（Y）;

toc;

实验结果（第一行为线性卷积，第二行为FFT）

L=8

M=8

时间已过0.232076秒。

时间已过0.087391秒。

L=16

M=16

时间已过0.000332秒。

时间已过0.000098秒。

L=32

M=32

时间已过0.000150秒。

时间已过0.000051秒。

L=64

M=64

时间已过0.289400秒。

时间已过0.000758秒。

L=256

M=256

时间已过0.000238秒。

时间已过0.000169秒。

L=512

M=512

时间已过0.000455秒。

时间已过0.000235秒。

L=1024

M=1024

时间已过0.001132秒。

时间已过0.000922秒。

分析：

可见在在相同长度下，快速卷积比线性卷积快差不多一倍的时间。

第二题

实验代码

实验结果：

（第一行为线性卷积，第二行为FFT）

L=2048

M=256

时间已过0.012061秒。

时间已过0.000581秒。

分析：

快速卷积比线性卷积快。

第三题

Lx=input（'L='）;

M=input（'M='）;

N=M+1;

x=1;

y=1;

fori=1:

L-1

x=[1x];

y=[ycos（0.2*pi*i）];

end

x=[x,zeros（1,N-1）];

t=zeros（1,M-1）;

Y=zeros（1,Lx+M-1）;

a=floor（Lx/N）;

tic

fork=0:

a

A=x（k*N+1:

k*N+N）;

y1=fft（A,Lx+M-1）;

y2=fft（h,Lx+M-1）;

y3=y1.*y2;

q=ifft（y3,Lx+M-1）;

Y（k*N+1:

k*N+M-1）=q（1:

M-1）+t（1:

M-1）;

Y（k*N+M:

k*N+N）=q（M:

N）;

t（1:

M-1）=q（N+1:

N+M-1）;

end

toc

实验结果：

L=2048

M=256

时间已过0.007656秒。

分析：

较第2题结果快。

第四题

实验代码：

Lx=input（'L='）;

M=input（'M='）;

N=M+1;

x=1;

y=1;

fori=1:

L-1

x=[1x];

y=[ycos（0.2*pi*i）];

end

M1=M-1;

L=N+M1;

h=[h,zeros（1,N-M）];

x=[zeros（1,M1）,x,zeros（1,N-1）];

a=floor（（Lx+M1-1）/（L））+1;

Y=zeros（1,N）;

tic

fork=0:

a-1

xk=x（k*L+1:

k*L+N）;

b=fft（xk,N）;

C=fft（h,N）;

Z=b.*C;

Y（k+1,:

）=ifft（Z,N）;

end

toc

L=2048

M=256

实验结果：

L=2048

M=256

时间已过0.002848秒。

分析：

较第2题结果快。

五、心得与体会

本次实验要求我们掌握利用FFT计算线性卷积的原理及具体实现方法，通过实验加深理解重叠相加法和重叠保留法并考察利用FFT计算线性卷积各种方法的适用范围。

本次实验让我切实看到了FFT算法的高效性及对于庞大的数据的处理实时性，同时对重叠保留法、重叠相加法也有了更深的认识，对以后的实验有了很好的理论基础，受益颇多。

实验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数字滤波器的优缺点及适用范围。

四、实验代码及实验结果

实验代码：

巴特沃斯：

fp=1000;

fs=1500;

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））;

Rp=1;

As=15;

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*N2）））;

[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;

xlabel（'\Omega（\pi）'）;ylabel（'|H（j\Omega）|'）;title（'°ÍÌØÎÖ˹Âö³åÂö³åÏìÓ¦²»±ä'）;

subplot（223）;

plot（w/pi,angle（H））;

gridon;

xlabel（'\Omega（\pi）'）;ylabel（'PhaseofH（j\Omega）|'）;title（'°ÍÌØÎÖ˹Âö³åÂö³åÏìÓ¦²»±ä'）;

[b2,a2]=butter（N2,OmegaC2,'s'）;

[H,w]=freqs（b2,a2）;

subplot（222）;

plot（w/pi,abs（H））;

axistight;

gridon;

xlabel（'\Omega（\pi）'）;ylabel（'|H（j\Omega）|'）;title（'°ÍÌØÎÖ˹˫ÏßÐԱ任'）;

subplot（224）;

plot（w/pi,angle（H））;

gridon;

xlabel（'\Omega（\pi）'）;ylabel（'PhaseofH（j\Omega）|'）;title（'°ÍÌØÎÖ˹˫ÏßÐԱ任'）;

实验结果：

切比雪夫Ⅰ型：

fp=1000;

fs=1500;

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））;

Rp=1;

As=15;

e=（10^（Rp/10）-1）^0.5;

N1=ceil（acos（（（As*As）-1）^0.5/e）/acos（Ws/Wp））;

N2=ceil（acos（（（As*As）-1）^0.5/e）/acos（ws/wp））;

[b1,a1]=cheby1（N1,Rp,Ws,'s'）;

[H,w]=freqs（b1,a1）;

subplot（221）;

plot（w/pi,abs（H））;

gridon;

xlabel（'\Omega（\pi）'）;ylabel（'|H（j\Omega）|'）;title（'ÇбÈÑ©·ò1Âö³åÂö³åÏìÓ¦²»±ä'）;

subplot（223）;

plot（w/pi,angle（H））;