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

1、 本科实验报告实验名称： 数字信号处理实验 课程名称：数字信号处理实验实验时间：任课教师：实验地点：4423实验教师：实验类型： 原理验证 综合设计 自主创新学生姓名：学号/班级：组 号：学 院：信息与电子学院同组搭档：专 业：信息工程成 绩：实验1 利用DFT分析信号频谱一、实验目的1.加深对DFT原理的理解。2.应用DFT分析信号频谱。3.深刻理解利用DFT分析信号频谱的原理，分析现实过程现象及解决办法。二、实验原理1、DFT和DTFT的关系 有限长序列的离散时间傅里叶变换在频率区间的N个等分点上的N个取样值可以由下式表示： 由上式可知，序列的N点DFT，实际上就是序列的DTFT在N个等间

2、隔频率点上样本。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分析连续时间信号的频谱 采用计算机分析连续时间信号的频谱，第一步就是把连续时间信号离散

3、化，这里需要进行连个操作：一是采样，二是截断。对于连续非周期信号，按采样间隔T进行采样，截取长度为M，那么对进行N点的频率采样，得到因此，可以将利用DFT分析连续非周期信号频谱的步骤归纳如下： （1）确定时域采样间隔T，得到离散序列； （2）确定截取长度M，得到M点离散序列，这里的为窗函数。 （3）确定频域采样点数N，要求。 （4）利用FFT计算离散序列的N点DFT，得到。 （5）根据式（2-6）由计算采样点的近似值。 采用上诉方法计算的频谱，需要注意如下三点问题：（1）频谱混叠。如果不满足采样定理的条件，频谱会很出现混叠误差。对于频谱无限宽的信号，应考虑覆盖大部分主要频率的范围。（2）栅栏效

4、应和频谱分辨率。使用DFT计算频谱，得到的结果只是N个频谱样本值，样本值之间的频谱是未知的，就像通过一个栅栏观察频谱，称为“栅栏效应”。频谱分辨率与记录长度成正比，提高频谱分辨率，就要增加记录时间。（3）频谱泄露。对于信号截断会把窗函数的频谱会引入到信号频谱中，造成频谱泄露。解决这问题的主要办法是采用旁瓣小的窗函数，频谱泄露和窗函数均会引起误差。 因此，要合理选取采样间隔和截取长度，必要时还需考虑适当的窗。 对于连续周期信号，我们在采用计算机进行计算时，也总是要进行截断，序列总是有限长的，仍然可以采用上诉方法近似计算。4、可能用到MATLAB函数与代码实验中的DFT运算可以采用MATLAB中提

5、供的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做

6、频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。4.利用DFT近似分析连续时间信号xt=e-0.1tu(t)的频谱（幅度谱）。分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定合适的参数。四、实验代码及实验结果实验1.1.1实验结果：实验代码： n=0:3; x=2 -1 1 1; w=-pi:0.01*pi:pi; X=x*exp(-j*n*w); subplot(211); plot(w,abs(X); title(幅度);xlabel(w);ylabel(|X|); axis tight; subplot(212); plot(w,angle(X); title(相位

7、);xlabel(w);ylabel(Angle（X）); axis tight;实验1.1.2实验结果：实验代码： n=0:3;x=2 -1 1 1;w=-pi:0.01*pi:pi;X=x*exp(-j*n*w);subplot(211);plot(w,abs(X);title(幅度);xlabel(w);ylabel(|X|);axis tight;hold on;subplot(212);plot(w,angle(X);title(相位);xlabel(w);ylabel(Angle（X）);axis tight;hold on;H=fft(x);subplot(211); stem(

8、n,abs(H),filled); subplot(212); stem(n,angle(H),filled);实验1.1.3实验结果：实验代码： x=2 -1 1 1 zeros(1,60); X=fft(x); subplot(211); n=0:63; stem(n,abs(X),filled); title(幅度);xlabel(n);ylabel(|X|); subplot(212); stem(n,angle(X),filled);title(相位);xlabel(n);ylabel(angle(X);分析：可以由DFT计算DTFT。通过补零加长序列，提高采样密度，可以由DFT近似

9、计算DTFT。实验1.2.1实验结果：实验代码：n=0:10;x=cos(0.48*n*pi)+cos(0.52*n*pi);X=fft(x);subplot(211);stem(n,abs(X),filled);title(幅度);xlabel(n);ylabel(|X|);subplot(212);stem(n,angle(X),filled);title(相位);xlabel(n);ylabel(angle(X);补零加长：实验结果：实验代码：h=x zeros(1,89);H=fft(h);n=0:99;subplot(211);stem(n,abs(H),filled);title(

10、幅度);xlabel(n);ylabel(|H|);subplot(212);stem(n,angle(H),filled);title(相位);xlabel(n);ylabel(angle(H);实验1.2.2实验结果：实验代码：n=0:100; x=cos(0.48*n*pi)+cos(0.52*n*pi); X=fft(x); subplot(211);stem(n,abs(X),filled);title(幅度);xlabel(n);ylabel(|X|);subplot(212);stem(n,angle(X),filled);title(相位);xlabel(n);ylabel(a

11、ngle(X);分析：可以通过增大截取长度和增加补零的个数来提高频谱分辨率，但是补零不能够增加分辨力。实验1.3实验结果：实验代码：f1=1;f2=2;f3=3;Fs=500;Tp=1;t=0:1/Fs:Tp;x=0.15*sin(2*pi*f1*t)+sin(2*pi*f2*t)-0.1*sin(2*pi*f3*t);N=501;F=Fs/N;f=0:Fs/(N-1):Fs;X=fft(x,N)；stem(f,abs(X)/250,filled);axis(0,4,0,1);xlabel(f/Hz);title(Magnitude);分析：通过选取合适的采样周期，可以完整的恢复出原信号的频谱

12、波形。实验1.4Tp=1s，0n100实验结果：实验代码：n=0:1:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),filled);Tp=5s,0n100实验结果：实验代码：n=0:5:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),filled);Tp=25s，0n100实验结果：实验代码：n=0:25:100;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),filled);Tp=1s，0n50实验结果：实验代码：n=0:1:50;x=exp(-0.1*n);X=fft(x);stem(n,abs

13、(X),filled);Tp=1,0n25实验结果：实验代码：n=0:1:25;x=exp(-0.1*n);X=fft(x);stem(n,abs(X),filled);分析： 最终确定参数：tp=5s，0n100。五、心得与体会通过本次实验，我们掌握并加深对DFT原理的理解并且学会应用DFT分析信号频谱，在此基础上利用DFT分析信号频谱的原理，掌握了利用matlab分析现实问题的步骤及办法。并且，通过这次的实验，对信号序列有了更加深刻的认识，单纯的一个信号序列是没有意义的，只有配合他本身的时间序列才是一个完整的信号序列，才可以对其进行分析。实验2 利用FFT计算线性卷积（选作）一、实验目的1

14、.掌握利用FFT计算线性卷积的原理及具体实现方法。2.加深理解重叠相加法和重叠保留法。3.考察利用FFT计算线性卷积各种方法的适用范围。二、实验原理1.线性卷积与圆周卷积设x(n)为L点序列，h(n)为M点序列，x(n)和h(n)的线性卷积为 （3-1）的长度为L+M-1x(n)和h(n)的圆周卷积为 （3-2）圆周卷积与线性卷积相等而不产生交叠的必要条件为NL+M+1 (3-3)圆周卷积定理：根据DFT性质，x(n)和h(n)的N点圆周卷积的DFT等于它们的DFT的乘积： (3-4)2.快速卷积快速卷积发运用圆周卷积实现线性卷积，根据圆周卷积定理利用FFT算法实现圆周卷积。可将快速卷积运算的

15、步骤归纳如下：(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

16、(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)。重叠相加

17、法：设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),0nL和

18、h(n)=cos(0.2n),0nM的线性卷积完成以下实验内容。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时计算线性卷积

19、的时间,与2题的结果进行比较。四、实验代码及实验结果实验2.1实验代码：for M=8 16 32 64 256 512 1024;L=M;n=0:1:L-1;x=ones(1,L);h=cos(0.2.*pi.*n);disp(L=M=);disp(M)disp(线性卷积用时：)ticy1=conv(x,h);tocdisp(快速卷积用时：);ticX1=fft(x);H=fft(h);Y=X1.*H;y=ifft(Y);tocend实验结果：L=M= 8线性卷积用时：时间已过 0.000122 秒。快速卷积用时：时间已过 0.000042 秒。L=M= 16线性卷积用时：时间已过 0.00

20、0045 秒。快速卷积用时：时间已过 0.000026 秒。L=M= 32线性卷积用时：时间已过 0.000051 秒。快速卷积用时：时间已过 0.000030 秒。L=M= 64线性卷积用时：时间已过 0.000050 秒。快速卷积用时：时间已过 0.000030 秒。L=M= 256线性卷积用时：时间已过 0.000089 秒。快速卷积用时：时间已过 0.000066 秒。L=M= 512线性卷积用时：时间已过 0.000540 秒。快速卷积用时：时间已过 0.001388 秒。L=M= 1024线性卷积用时：时间已过 0.000617 秒。快速卷积用时：时间已过 0.000115 秒。分

21、析：可见在在相同长度下，快速卷积比线性卷积快差不多一倍的时间。实验2.2实验代码：L=2048;M=256;n1=0:1:L;n2=0:1:M;x=ones(1,L);h=cos(0.2.*pi.*n2);disp(线性卷积：)ticy1=conv(x,h);tocdisp(快速卷积：)N=M+L-1;ticX1=fft(x,N);H=fft(h,N);Y=X1.*H;y=ifft(Y,N);toc实验结果： L=2048，M=256线性卷积：时间已过 0.000437 秒。快速卷积：时间已过 0.000407 秒。L=4096，M=256线性卷积：时间已过 0.000445 秒。快速卷积：时

22、间已过 0.001395 秒。分析：当序列长度较短时，快速卷积比线性卷积快；当序列长度过长时，快速卷积比线性卷积慢了许多。实验2.3实验代码：L=2048;M=256;disp(L=);disp(L);disp(M=);disp(M);x=ones(1,L);n=0:1:M-1;h=cos(0.2.*pi.*n);Li=M;Num=L/Li;Ni=Li+M-1;N_yu=L+M-1;y = zeros(1,N_yu);disp(重叠相加法卷积：);ticH = fft(h,Ni);for i=1:Num;i_low=(i-1)*M+1;xi=x(i_low:i_low+M-1);Xi=fft(

23、xi,Ni);Yi=Xi.*H;yi=ifft(Yi,Ni);y(i_low:i_low+Ni-1)=y(i_low:i_low+Ni-1)+yi;endtoc实验结果：L= 2048M= 256重叠相加法卷积时间时间已过 0.000626 秒。分析：较第2题结果快。实验2.4实验代码：L=2048;M=256;disp(L=);disp(L);disp(M=);disp(M);x=ones(1,L);n=0:1:M-1;h=cos(0.2.*pi.*n);Li=M;Num=L/Li;Ni=Li+M-1;N_yu=L+M-1;y=zeros(1,N_yu);x_add=zeros(1,M-1)

24、,x;disp();ticH = fft(h,Ni);for i=1:Num;i_low=(i-1)*M+1;xi=x_add(i_low:i_low+M-1);Xi=fft(xi,Ni);Yi=Xi.*H;yi=ifft(Yi,Ni);y(i_low:i_low+M-1)=yi(M:Ni);endxi=x_add(i*M+1):N_yu);Xi=fft(xi,Ni);Y=Xi.*H;yi=ifft(Y,Ni);y(i*M+1):N_yu)=yi(M:Ni-1);toc实验结果：L= 2048M= 256重叠保留法卷积时间时间已过 0.000508 秒。分析：较第2题结果快。五、心得与体会本次

25、实验要求我们掌握利用FFT计算线性卷积的原理及具体实现方法，以及通过对快速卷积和线性卷积运算速度的比较，更加直观的看到利用FFT快速卷积的优点。本次实验让我切实看到了FFT算法的高效性及对于庞大的数据的处理实时性，同时对重叠保留法、重叠相加法也有了更深的认识，对以后的实验有了很好的理论基础，受益颇多。实验3 IIR数字滤波器设计一、实验目的1.掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法。2.加深理解数字滤波器和模拟滤波器之间的技术指标转化。3.掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围。二、实验内容1.设采样频率为fs =4kHz，采用

26、脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3dB截止频率为fc=1kHz。2.设采样频率为fs=10kHz，设计数字低通滤波器，满足如下指标 通带截止频率：fp=1kHz，通带波动：Rp=1dB 阻带截止频率：fst=1.5kHz，阻带衰减：As=15dB要求分别采用巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。结合实验结果，分别讨论采用上述设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。四、实验代码及实验结果巴特沃斯：实验代码：fs=10*1000;fp=1*1000;fst=1.5*1000

27、;Rp=1;As=15;Wp = 2*pi*fp/fs;Ws = 2*pi*fst/fs;N = ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws);Wc = Wp /(10(Rp/10)-1)(1/(2*N);b,a = butter(N,Wc/pi,s);bz,az = impinvar(b,a,1)w=0:500*pi/500;H,w=freqz(bz,az); %频率响应subplot(221);plot(w/pi,abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|);subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H);grid on;xlabel(omega(pi);ylabel(|H(ejomega)|(dB);subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel(omega(pi);ylabel(Phase of H(ejomega)(pi);grd=grpdelay(bz,az,w); %群延时响应subplot(224);plot(w/pi,grd);grid on;xlabel(omega(pi);ylabel(Gr