数字信号处理课程研究性学习报告.doc

数字信号处理课程研究性学习报告姓名 学号 同组成员 指导教师 申艳 时间 2015年6月10日星期三 题目一：基本概念和技能学习报告【目的】(1) 掌握离散信号和系统时域、频域和z域分析中的基本方法和概念；(2) 学会用计算机进行离散信号和系统时域、频域和z域分析。(3) 培养学生自主学习能力，以及发现问题、分析问题和解决问题的能力。M1-1已知，以抽样频率对上述三个信号进行抽样，在同一张图上画出，和及其抽样点，对所得结果进行讨论。【题目目的】1. 掌握抽样的基本概念；2. 学会MATLAB中对信号抽样的方法。【仿真结果】【结果分析】通过仿真抽样前后的信号图像我们可知，g1=cos(6*pi*t),g2=cos(14*pi*t),g3=cos(26*pi*t)虽然它们信号时域表达式不同，但在以抽样频率fm=10Hz抽样后，各离散信号具有相同的时域表达式，所以当抽样频率过小时，所得离散信号可能存在较大失真，无法反映原信号的特征，所以在给定一个信号后，我们应给出合适的抽样频率，抽样频率过小时失真较多，抽样频率过大时又造成了浪费，抽样信号符合我们给定的要求即可，这样得到的信号才能在最大程度保留原信号的特征的基础上节约资源。【仿真程序】k1=0:0.01:2;k2=0:0.1:2;g1=cos(6*pi*k1);g2=cos(6*pi*k2);subplot(3,2,1);plot(k1,g1);axis(0,2,-1.5,1.5);title(g1=cos(6*pi*t) 原信号);xlabel(t);ylabel(g1(t);subplot(3,2,2);stem(k2,g2);axis(0,2,-1.5,1.5);title(g1=cos(6*pi*t) 抽样频率fm=10Hz);xlabel(k);ylabel(g1k); k3=0:0.01:2;k4=0:0.1:2;g3=cos(14*pi*k3);g4=cos(14*pi*k4);subplot(3,2,3);plot(k3,g3);axis(0,2,-1.5,1.5);title(g2=cos(14*pi*t) 原信号);xlabel(t);ylabel(g2(t);subplot(3,2,4);stem(k4,g4);axis(0,2,-1.5,1.5);title(g2=cos(14*pi*t) 抽样频率fm=10Hz);xlabel(k);ylabel(g2k); k5=0:0.001:2;k6=0:0.1:2;g5=cos(26*pi*k5);g6=cos(26*pi*k6);subplot(3,2,5);plot(k5,g5);axis(0,2,-1.5,1.5);title(g3=cos(26*pi*t) 原信号);xlabel(t);ylabel(g3(t);subplot(3,2,6);stem(k6,g6);axis(0,2,-1.5,1.5);title(g3=cos(26*pi*t) 抽样频率fm=10Hz);xlabel(k);ylabel(g3k);M1-2利用MATLAB的filter函数，求出下列系统的单位脉冲响应，并判断系统是否稳定。讨论本题所获得的结果。【题目目的】1. 掌握LTI系统单位脉冲响应的基本概念、系统稳定性与单位脉冲响应的关系；2. 学会filter函数的使用方法及用filter函数计算系统单位脉冲响应；3. 体验有限字长对系统特性的影响。【仿真结果】下面是两个系统的单位脉冲响应：下面是两个系统的零极点分布图：极点1 0.9430 0.9020极点2 1.0000 0.8500【结果分析】根据以上仿真结果可知，两个系统函数H1，H2虽然只是分母系数差了一点，但最终实现时系统的单位脉冲响应会出现很大的不同，系统的极点也发生了不小的变化，这就造成了系统的失真。比如，完成数字滤波器设计后，具体实现时，可能会出现较大的误差，甚至使设计的稳定系统变成不稳定系统，从而无法达到设计的要求。引起这些误差的根本原因在于寄存器（存储单元）的字长有限。误差的特性与系统的类型、结构形式、数字的表示法、运算方式及字的长短有关。在通用计算机上，字长较长，量化步很小，量化误差不大。但在专用硬件，如FPGA，实现数字系统时，其字长较短，就必须考虑有限字长效应了。【问题探究】已知LTI系统的系统函数，有哪些计算系统单位脉冲响应方法，比较这些方法的优缺点。 filter函数，可计算出差分方程的零状态响应，既可以用来求yk,也可以求出hk，是最基本的方法，也是最复杂的方法;impulse函数，只是用来实现冲击响应的，是最方便的方法；conv函数，是用来计算卷积的，可以用来求yk，这种方法也较麻烦，间接求出单位脉冲响应，不能直接求出单位脉冲响应。【仿真程序】b1=1;a1=1,-1.845,0.850586;b2=1;a2=1,-1.85,0.85;x=1,zeros(1,49);y1=filter(b1,a1,x);subplot(2,1,1);stem(y1);xlabel(x);ylabel(h1k);axis(0 50 0 8);y2=filter(b2,a2,x);subplot(2,1,2);stem(y2);xlabel(x);ylabel(h2k);subplot(2,1,1);r1,p1,m1=residuez(b1,a1); disp(极点1);disp(p1);zplane(b1,a1);title(H1);subplot(2,1,2);r2,p2,m2=residuez(b2,a2); disp(极点2);disp(p2);zplane(b2,a2);title(H2);M1-3(1)利用MATLAB语句x=firls(511,0 0.4 0.404 1,1 1 0 0）产生一个长度为512的序列xk，并画出该序列的幅度频谱。(2) 已知序列，分别画出时序列yk的幅度频谱。解释所得到的结果。【题目目的】1. 学会用MATLAB函数freqz计算序列频谱；2. 掌握序列频谱的基本特性及分析方法。【温磬提示】只需知道MATLAB语句x=firls(511,0 0.4 0.404 1,1 1 0 0）产生一个长度为512的序列xk，该序列满足不需知道其他细节。用函数freqz计算该序列的频谱，在画幅度频谱时，建议用归一化频率。【仿真结果】改进后：【结果分析】yk=对信号进行Fourier变换：所以，yk的幅度频谱是xk幅度频谱左移与右移之和的一半，在平移的过程中，会出现混叠现象，然而这样的理论分析与实际的结果不同，比如w0=0.8pi时，本来值较高的地方实际的值却很低，原因如下：在xk变化较剧烈的点的值处与cos函数值为零的值相乘后消去，所以将高频分量丢失，所以频谱会出现理论与实际结果不同的现象。所以，是512的偶数倍的点均会出现这种现象，所以可以选择将xk变成3000个点的序列，即可解决这种现象。【问题探究】有部分的计算结果可能与理论分析的结果不一致，分析出现该现象的原因，给出解决问题方法并进行仿真实验。【仿真程序】x=firls(511,0 0.4 0.404 1,1 1 0 0);b=x;a=1;w=linspace(-pi,pi,512);y=freqz(b,a,w);subplot(3,2,1:2);plot(w/pi,abs(y);axis(-2,2,-1,2);title(生成序列xk的幅度频谱);k=0:511;y1=x.*cos(0.4*pi*k);g1=freqz(y1,a,w);subplot(3,2,3);plot(w/pi,abs(g1);xlabel(w0=0.4pi);ylabel(yk=xkcos(w0k)的幅度频谱);y2=x.*cos(0.8*pi*k);g2=freqz(y2,a,w);subplot(3,2,4);plot(w/pi,abs(g2);xlabel(w0=0.8pi);ylabel(yk=xkcos(w0k)的幅度频谱);y3=x.*cos(0.9*pi*k);g3=freqz(y3,a,w);subplot(3,2,5);plot(w/pi,abs(g3);xlabel(w0=0.9pi);ylabel(yk=xkcos(w0k)的幅度频谱);y4=x.*cos(pi*k);g4=freqz(y4,a,w);subplot(3,2,6);plot(w/pi,abs(g4);xlabel(w0=pi);ylabel(yk=xkcos(w0k)的幅度频谱);改进后的程序：x=firls(2999,0 0.4 0.404 1,1 1 0 0);b=x;a=1;w=linspace(-pi,pi,512);y=freqz(b,a,w);subplot(3,2,1);plot(w/pi,abs(y);axis(-2,2,-1,2);title(生成序列xk的幅度频谱);k=0:2999;y1=x.*cos(0.4*pi*k);subplot(3,2,2);plot(k,x);title(xk的时域图像);g1=freqz(y1,a,w);subplot(3,2,3);plot(w/pi,abs(g1);xlabel(w0=0.4pi);ylabel(yk=xkcos(w0k)的幅度频谱);y2=x.*cos(0.8*pi*k);g2=freqz(y2,a,w);subplot(3,2,4);plot(w/pi,abs(g2);xlabel(w0=0.8pi);ylabel(yk=xkcos(w0k)的幅度频谱);y3=x.*cos(0.9*pi*k);g3=freqz(y3,a,w);subplot(3,2,5);plot(w/pi,abs(g3);xlabel(w0=0.9pi);ylabel(yk=xkcos(w0k)的幅度频谱);y4=x.*cos(pi*k);g4=freqz(y4,a,w);subplot(3,2,6);plot(w/pi,abs(g4);xlabel(w0=pi);ylabel(yk=xkcos(w0k)的幅度频谱);M1-4已知，当，比较的幅度响应。从中你能得出什么结论？ 【题目目的】1. 掌握低通滤波器、高通滤波器、FIR滤波器和IIR滤波器的概念；2. 学会使用freqz函数。【仿真结果】【结果分析】由H1k的幅度响应图线可知，H1代表的是数字低通FIR滤波器，当=0.8时，H2代表的是数字低通IIR滤波器，当=-0.8时，H2代表的系统是数字IIR滤波器，并且高频分量幅度响应很大。当=0.8时，H3代表的是数字低通IIR滤波器，是H1和H2级联的结果，当=-0.8时，H3代表的是数字全通IIR滤波器，幅度响应近似为1，是H1和H2级联的结果。由此可知，当参数不同时，系统会有较大的不同。【仿真程序】b1=1,1;a1=2;w=linspace(-pi,pi,512);y1=freqz(b1,a1,w);subplot(3,1,1);plot(w/pi,abs(y1);axis(-1.5,1.5,-0.5,1.5);title(H1(z)的幅度响应);xlabel(归一化频率w/pi);ylabel(|H1|);subplot(3,1,2);b21=0.2;a21=1 -0.8;b22=1.8;a22=1 0.8;y21=freqz(b21,a21,w);y22=freqz(b22,a22,w);plot(w/pi,abs(y21),r,w/pi,abs(y22);title(H2(z)的幅度响应);xlabel(归一化频率w/pi);ylabel(|H2|);axis(-1.5,1.5,-2,12);b31=conv(b1,b21);a31=conv(a1,a21);b32=conv(b1,b22);a32=conv(a1,a22);y31=freqz(b31,a31,w);y32=freqz(b32,a32,w);subplot(3,1,3);plot(w/pi,abs(y31),r,w/pi,abs(y32);title(H3(z)的幅度响应);xlabel(归一化频率w/pi);ylabel(|H3|);axis(-1.5,1.5,-0.5,1.5);M1-6已知一因果系统的H(z) 为试求出和H(z)具有相同幅度响应的最小相位系统Hmin(z)和最大相位系统Hmax(z)。画出并比较H(z)、Hmin(z) 和Hmax(z)的相位响应。【题目目的】1. 掌握全通滤波器的基本特征和特性；2. 学会计算具有相同幅度响应的最小相位系统Hmin(z)和最大相位系统Hmax(z)。【温磬提示】在画系统的相位响应时，要把三个系统的相位响应画在一张图上，比较其解卷绕(unwrap)后的相位特点。【仿真结果】手算结果：Hmin(z)=，Hmax(z)=计算机自动求解的结果：【结果分析】根据图像可以看出，在具有相同幅度响应的因果稳定系统中，最小相位系统具有最具有最小的相位滞后，最大相位系统具有最大的相位滞后，原系统的相位滞后介于两者之间。【问题探究】手算时，如何找出最小相位系统Hmin(z)和最大相位系统Hmax(z)？用计算机自动求解时，如何找出最小相位系统Hmin(z)和最大相位系统Hmax(z)？ 你所用的算法是一样的吗？ 原理类似。【仿真程序】手算后：b1=1.1,1.99,0.9;a1=1,1.55,0.6;b2=0.9,1.99,1.1;a2=a1;b3=1,2,0.99;a3=a2;w=linspace(0,pi,1000);h1=freqz(b1,a1,w);h2=freqz(b2,a2,w);h3=freqz(b3,a3,w);plot(w,angle(h1),w,unwrap(angle(h2),r,w,angle(h3),g);axis(0,pi,-8,2);xlabel(w);ylabel(相位响应);计算机自动求解：b1=1,2,0.99;a1=1,1.55,0.6;w=linspace(0,pi,1000);h1=freqz(b1,a1,w);j1=roots(b1);j2=j1;j3=j1;if j1(1)1|j1(1)1|j1(2)N时则将原序列截短为N点序列再进行计算，若NM则将序列补零后在做N点的DFT。【仿真结果】L=32L=64L=128L=256L=512L=32L=128【结果分析】对序列补零后再做DFT相当于增加了DFT的点数，即频域抽样的点数，而原离散序列没有改变，其傅里叶变换结果也没有改变，同时若DFT点数太少则获得的频谱信息过少，有可能会使得重要的频率信息丢失。由结果可知，DFT点数越多，产生的离散谱中含有的信息也就越多，得到的频谱能更好的反应原连续谱中的信息。在对离散序列用DFT做谱分析时，应当适当增加DFT的点数，以减小栅栏效应。DFT点数越多，则L越大，即fsam/N越小，那么显示分辨率就越高。 【自主学习内容】已知幅度谱，用matlab求谱峰所对应的频率值。利用：a,b=max(X);其中X是一个向量，返回的b值对应的是向量X中的最大值，返回的a值为最大值b对应的索引。【阅读文献】1陈后金，薛健，胡健. 数字信号处理M. 北京：高等教育出版社，2006 .【发现问题】 (专题研讨或相关知识点学习中发现的问题)：随着DFT点数的增加，所获得的频谱图像与实际值误差减小，而谱峰值更精确，且其对应的谱峰频率变大。而在用matlab模拟近似频谱时，由于是离散的点构成的，无法正常表现出真实频谱的图像【问题探究】【仿真程序】 k=0:31; x=sin(0.2*pi*k); L=0:31; X=fft(x,32);%对函数进行32点的DFT抽样，计算频谱 plot(2*L/32,abs(X)%abs()是求绝对值 title(L=32时序列的频谱) grid on %加网格线 Xmax=max(abs(X);%找到谱峰 y1=find(abs(X)=Xmax)-1;%找到谱峰对应的点,因为点从0开始所以要减1 m1=y1*2/32; disp(Xmax);14.5561 disp(m1);0.1875 1.8125k=0:31;x=sin(0.2*pi*k);L=0:63;X=fft(x,64);plot(2*L/64,abs(X)title(L=64时序列的频谱)grid on Xmax=max(abs(X);%找到谱峰 y1=find(abs(X)=Xmax)-1;%找到谱峰对应的点,因为点从0开始所以要减一 m1=y1*2/64; disp(Xmax); 14.5561 disp(m1); 0.0938 0.9063 k=0:31;x=sin(0.2*pi*k);L=0:127;X=fft(x,128);plot(2*L/128,abs(X)title(L=128时序列的频谱)grid on Xmax=max(abs(X);%找到谱峰 y1=find(abs(X)=Xmax)-1;%找到谱峰对应的点,因为点从0开始所以要减一 m1=y1*2/128;disp(Xmax); 14.5561 disp(m1); 0.0469 0.4531 k=0:31;x=sin(0.2*pi*k);L=0:255;X=fft(x,256);plot(2*L/256,abs(X)title(L=256时序列的频谱)grid on Xmax=max(abs(X);%找到谱峰 y1=find(abs(X)=Xmax)-1;%找到谱峰对应的点,因为点从0开始所以要减一 m1=y1*2/256;disp(Xmax); 14.5561 disp(m1); 0.0234 0.2266 k=0:31;x=sin(0.2*pi*k);L=0:511;X=fft(x,512);plot(2*L/512,abs(X)title(L=512时序列的频谱)grid on Xmax=max(abs(X);%找到谱峰 y1=find(abs(X)=Xmax)-1;%找到谱峰对应的点,因为点从0开始所以要减一 m1=y1*2/512;disp(Xmax); 14.5561 disp(m1); 0.0117 0.1133M2-3 已知一离散序列为 x k=AcosW0k+Bcos ( (W0+DW)k)。用长度N=64的哈明窗对信号截短后近似计算其频谱。试用不同的A和B的取值，确定用哈明窗能分辨的最小的谱峰间隔中c的值。【题目分析】本题讨论用哈明窗计算序列频谱时的频率分辨率。Hamming窗函数的幅值有中心向两端逐渐减弱，因而其高频分量明显减小，频谱中旁瓣的幅度较小，主瓣峰值与第一个旁瓣峰值相对衰减很大，hamming窗以增加主瓣宽度来降低旁瓣能量，用hamming窗极端频谱时要求能分辨的谱峰的间隔c/Tp=c*fs/N。针对不同的AB值，调整频率的差值来比较哈明窗所分辨的最小的谱峰。利用控制变量的思想，首先固定AB值，比较不同频率所造成的影响。然后固定频率，比较不同的AB值所造成的影响。解题思路：XK加窗截断，得到XnK，对其进行DFT计算，然后画出频谱。【仿真结果】Dw1=pi/32,C=1Dw2=pi/16,C=2Dw3=3*pi/32,C=3 dw4=pi/8,C=4 谱峰间隔相等时AB相差较大与AB近似相等两种情况：A：B=1:10A:B=1：1【结果分析】将实验结果与教材中定义的窗函数的有效宽度相比较，发表你的看法。（1） 当AB近似时，N=3时就能有效的区别能分辨出谱峰间隔。C=2是教材中定义的哈明窗的有效宽度，在实际中有时候是区别分开的。从图像可以看出C=3的时候就能完全分辨出来，C值越大，显示的额频谱越接近真实值。C=4是哈明窗的主瓣宽度，就能完全区别谱峰间隔。（2） 当AB相差较大的时候，采用相同的分析方法。我们可以确定大约C=3能够有效的分辨出谱峰间隔。【自主学习内容】几种常见的窗函数的特点与应用。有关DFT频谱分析，对连续离散信号的分析处理，频谱分辨率，谱峰间隔的知识。【阅读文献】1陈后金，薛健，胡健. 数字信号处理M. 北京：高等教育出版社，2006 .【问题探究】 在离散序列频谱计算中为何要用窗函数？用不同的窗函数对计算结果有何影响？与矩形窗相比哈明窗有何特点?如何选择窗函数？【仿真程序】A=B时，比较不同的谱峰间隔 N=64; k=0:63; L=0:511; A=1;B=1; dw1=pi/32; x1=A*cos(pi*k./4)+B*cos(pi*k./4+dw1*k);wh=(hamming(N); x1=x1.*wh; y1=fft(x1,512); figure;plot(L/512,abs(y1);xlabel(归一化频率); ylabel(幅值); title(dw1=pi/32,C=1) dw2=pi/16;x2=A*cos(pi*k./4)+B*cos(pi*k./4+dw2*k);wh=(hamming(N);x2=x2.*wh;y2=fft(x2,512);figure;plot(L/512,abs(y2); xlabel(归一化频率);ylabel(幅值);title(dw2=pi/16,C=2)dw3=3*pi/32;x3=A*cos(pi*k./4)+B*cos(pi*k./4+dw3*k);wh=(hamming(N);x3=x3.*wh;y3=fft(x3,512);figure;plot(L/512,abs(y3); xlabel(归一化频率);ylabel(幅值);title(dw3=3*pi/32,C=3) dw4=pi/8;x4=A*cos(pi*k./4)+B*cos(pi*k./4+dw4*k);wh=(hamming(N);x4=x4.*wh;y4=fft(x4,512);figure;plot(L/512,abs(y4); xlabel(归一化频率);ylabel(幅值);title(dw4=pi/8,C=4)A与B相差较大时 N=64;k=0:63;L=0:511;A=1;B=10;dw1=pi/32;x1=A*cos(pi*k./4)+B*cos(pi*k./4+dw1*k);wh=(hamming(N); x1=x1.*wh;y1=fft(x1,512);figure;plot(L/512,abs(y1);xlabel(归一化频率); ylabel(幅值);title(dw1=pi/32,C=1);M2-4 已知一离散序列为 xk=cos(W0k)+0.75cos(W1k), 0 k 63 其中W0=0.4p, W1=W0+p/64(1) 对xk做64点FFT, 画出此时信号的频谱。(2) 如果(1)中显示的谱不能分辨两个谱峰，是否可对(1)中的64点信号补零而分辨出两个谱峰。通过编程进行证实，并解释其原因 。(3) 给出一种能分辨出信号中两个谱峰的计算方案，并进行仿真实验。【题目分析】分析影响谱峰分辨率的主要因素，进一步认识补零在在频谱计算中的作用。分析影响谱峰分辨率的主要因数，进一步认识补领在频谱计算中的作用。影响谱峰分辨率的主要因数是序列的频率，即连续时间的抽样点数。【仿真结果】（1）（3） L=128L=512(3)64点信号改为128点信号，即能分辨出来。L=128L=256【结果分析】由结果可知，对信号增加零点不能提高频谱的分辨率；为了提高频谱的分辨率，可以增加序列的抽样点数，比如开始的时候，K的范围是0，64，后来K的范围是0，128，然后是0，256随着抽样点数的增加，频谱的分辨率逐渐提高。【自主学习内容】如何增加频谱分辨率【阅读文献】1陈后金，薛健，胡健. 数字信号处理M. 北京：高等教育出版社，2006 .【问题探究】离散信号的频谱具有周期性，会因抽样频率的不适而产生混叠，导致失真。【仿真程序】（1） k=0:63; L=64; f1=0.4*pi; f2=f1+pi/64; x=cos(f1*k)+0.75*cos(f2*k); x=x zeros(1,L-length(x); X=fftshift(fft(x); f=-pi:2*pi/L:pi-2*pi/L; plot(f/pi,abs(X) grid（2） k=0:63;L=128;f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fftshift(fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,abs(X)grid(3)k=0:127;L=128;f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fftshift(fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,abs(X); title(128点DFT);gridk=0:255;L=256;f1=0.4*pi;f2=f1+pi/64;x=cos(f1*k)+0.75*cos(f2*k);x=x zeros(1,L-length(x);X=fftshift(fft(x);f=-pi:2*pi/L:pi-2*pi/L;plot(f/pi,abs(X);title(256点DFT);gridM2-6 试用DFT近似计算高斯信号的频谱抽样值。高斯信号频谱的理论值为通过与理论值比较，讨论信号的时域截取长度和抽样频率对计算误差的影响。【题目分析】连续非周期信号频谱计算的基本方法。计算中出现误差的主要原因及减小误差的方法。利用连续非周期信号频谱与离散Fourier变换的关系。对于连续非周期信号，要对其进行频谱计算，需要先经过抽样，变为离散信号，再对离散信号的频谱在一个周期内抽样，即做DFT。由于题目所给信号时域无线长，所以为了在抽样后频谱不混叠，对信号进行截短即加抗混叠滤波器后在进行频谱分析，因为加滤波器后照成的误差远小于频谱混叠照成的误差。在对信号抽样时还需要考虑抽样频率，以使抽样后所得频谱不会发生混叠。本题以d=0.05为例，讨论fs、N变化带来的影响；再讨论d不同带来的影响。【仿真结果】fs=1,N=4,d=0.05fs=1,N=6,d=0.05fs=1,N=8,d=0.05fs变化：fs=0.1;N=128;d=0.05fs=0.2;N=128;d=0.05fs=0.4;N=128;d=0.05d变化时：fs=1;N=128;d=0.05fs=1;N=128;d=0.5fs=1;N=128;d=2【结果分析】由于信号及频谱都有理论表达式，在进行误差分析时希望给出一些定量的结果。当N逐渐变大时，谱线逐渐变为谱峰，而且谱峰的宽度随在一定范围内变化时逐渐变窄；当f变化时，对DFT分析信号频谱的精度影响较大，因为其直接影响频谱混叠的精度；当d逐渐变大时，N合适大小时，需要改变抽样频率来使的近似值逼近理论曲线。【自主学习内容】利用DFT分析连续非周期信号的频谱的相关概念及分析过程中出现的现象。【阅读文献】1陈后金，薛健，胡健. 数字信号处理M. 北京：高等教育出版社，2006 .【发现问题】 (专题研讨或相关知识点学习中发现的问题)：信号关于t=0偶对称，在计算过程中如何处理这个问题。【问题探究】N，fs，d的取值问题【仿真程序】 fs=1; N=4; T=1/fs;ws=2*pi*fs; d=0.05; L=512; g=exp(-d*(-N+1:N-1)*T.*(-N+1:N-1)*T); X=fftshift(fft(g,L); w=(-ws/2+(0:L-1)*ws/L)/(2*pi); wt=linspace(-2*pi,2*pi,2048); G=1/T*sqrt(pi/d).*exp(-wt.*wt)/(4*d); plot(w,abs(X),wt/(2*pi),abs(G); legend(近似值,理论值) grid扩展题M2-7 本题研究连续周期信号频谱的近似计算问题。 周期为T0的连续时间周期信号x(t)可用Fourier级数表示为其中X(nw0)称为连续时间周期信号x(t)的频谱函数。称为信号的基频（基波），称为信号的谐波。如果信号x(t)函数表达式已知，则可由积分得出信号的频谱。如果信号x(t)函数表达式未知，或者x(t)函数表达式非常复杂，则很难由积分得信号的频谱。本题的目的就是研究如何利用DFT近似计算连续时间周期信号的频谱。(1)若在信号x(t)的一个周期T0内抽样N个点，即， T为抽样周期(间隔)，可获得序列xk试分析序列xk的DFT与连续时间周期信号x(t)的频谱X(nw0)的关系；(2)由(1)的结论，给出由DFT近似计算周期信号频谱X(nw0)的方案；(3)周期信号x(t)的周期T0=1，x(t)在区间0，1的表达式为x(t)=20t2(1-t)4cos(12pt)(a)试画出信号x(t)在区间0，1的波形；(b)若要用10次以内的谐波近似表示x(t)，试给出计算方案，并计算出近似表示的误差。讨论出现误差的原因及减小误差的方法。 【理论推导】DFT计算所得结果Xm与连续周期信号频谱X(nw0)的关系。【计算方案】根据理论推导结果设计近似计算方案。分析产生误差的主要原因。Xm是连续周期信号频谱X(n0)的周期化，在乘上比例因子N所得，由于对于一般连续周期信号的频谱，当n时，X（n0）0,所以可以用Xm来近似较为低次的谐波分量，但仍会有一些高频谐波分量在周期化是发生混叠。如果信号是带限信号，但抽样时不满足抽样定理，也同样会产生混叠误差。用矩形窗截取了一个周期长度上的信号抽样值，因此会产生频率泄露。如果信号不是带限信号，在分析其频谱的时候，要截断成时限的谱，此时会产生截断误差。在对信号做DFT时，相当于在时域做了周期化，则对应的频谱进行离散化，就会产生栅栏效应。用有限次谐波合成信号，丢失的谐波有一定的误差在信号的重建是，谐波分量选取的次数多少也对误差有着较大的影响。对于一般的周期信号来说，通常是不带限的，而我们只能选择次数较低的谐波分量（占原有信号的主要能量）进行信号重建，这是不可避免的会发生泄漏，产生误差。【扩展分析】如果周期信号x(t)是带限信号，即信号的最高频率分量为Mw0(是正整数)，试确定在一个周期内的最少抽样点N，使得在频谱的计算过程当中不存在混叠误差。与抽样定理给出的结论比较，发表你的看法。由上面分析可知，频谱以N0周期化，原信号x(t)带限，最高频率分量为M0，当N2M时，频谱周期化不会发生混叠，可以看出这个关系式与抽样定理十分相像。周期信号同样遵循抽样定理，当抽样频率fsam2fm时，其中fsam=N/t0,fm=M0/2=M/T0，于是有N2M为不发生混叠的条件。【仿真结果】【结果分析】讨论DFT点数对近似计算的影响，讨论所取谐波项的多少对近似计算的影响。误差分析要给出定量的结果，如平均误差，最大误差等。【仿真程序】N=32; T0=1;n=10; T=1/N; ts=(0:N-1)*T; xs=20.*ts.*ts.*(1-ts).4).*cos(12*pi*ts); X=fft(xs,N);C=X/N;t=linspace(0,1,1024);xt=20.*t.*t.*(1-t).4).*cos(12*pi*t);figure(1);plot(t,xt);xlabel(t)ylabel(xt)title(xt的时域波形);grid onfigure(2)k=0:N-1;stem(k,abs(X);xlabel(k)ylabel(X)title(xs的N点DFT);grid onfigure(3);xe=C(1)*ones(1,length(t); for k=1:n-1 xe=xe+C(k+1)*exp(j*k*2*pi/