信号系统关于栅栏效应和吉布斯效应验证

1、佳木斯大学课程设计报告基于MATLAB的信号吉布斯效应和栅栏效应的验证学 院 信息电子技术专 业 班 级 学 籍 号 姓 名 指导教师 佳 木 斯 大 学2015年7月5日佳木斯大学课程设计报告第一章 MATLAB简介1.1 MATLAB简介MATLAB是美国Mathworks公司推出的数学工具软件， 它是一种直观、高效的计算机语言，同时又是一个科学计算平台。它为数据分析和数据可视化、算法提供了最核心的数学和高级图形工具。多达几百个数学函数和工程函数，极大地降低了对使用者的数学基础和计算机语言知识的要求，编程效率和计算效率极高，初学者可在很短的时间就能初步掌握它。目前，MATLAB 已经成为国

2、际上最流行的科学与工程计算软件工具。它集科学计算、图像处理、声音处理于一身，是一个高度的集成系统，有良好的用户界面，并有良好的帮助功能。MATLAB不仅流行于控制界，在机械工程、生物工程、语音处理、图像处理、信号分析、计算机技术等各行各业中都有极广泛的应用。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。它的MATLAB的基本数据单位是矩阵，

3、它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C、FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，因此使MATLAB成为一个强大的数学软件。MATLAB语言的特点：编程效率高、用户使用方便、扩充能力强 、语句简单、内涵丰富、高效方便的矩阵和数组运算、方便的绘图功能。第二章 离散傅立叶变换2.1 傅立叶变换FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱

4、提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。现在就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍。采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率

5、值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。傅里叶变换常用的两种傅里叶

6、变换是：周期序列的离散傅里叶级数(DFS)和非周期序列的傅里叶变换(DTFT)，其表示式分别为： (2-1) (2-2)当用数字计算机对信号进行频谱分析时，要求信号必须以有限长度的离散值作为输入，而计算所得的频谱值自然也是有限、离散的。但其时间变量和频率变量又同时具有周期性；DTFT变换满足时间自变量的有限长度，但其频率变量为连续形式。可见，这两种变换都难以实际应用。离散傅里叶变换(DFT)是对序列的傅里叶变换(DTFT)的频域采样，运用频域采样方法同样可以推导出DFT的变换公式。对于任意函数可分解为无穷多个不同频率正弦信号的和, 即频谱分析。 2.2 离散FT(DFT)离散傅里叶变换使数字计

7、算机和信号与系统分析紧密地结合,虽然离散傅里叶变换在信号处理中发挥重要作用,但它的定义并不严格,离散傅里叶变换的定义是从离散傅里叶级数引入的,而离散傅里叶级数是由混迭误差的概念推出并直接定义,任意长序列傅里叶变换至今还没有定义,对常遇到的非周期序列, 有限长或无限, 只能作DTFT, 即连续频谱X(), 模拟在计算机上做数值计算,实际中, 把N点序列视为一周期序列的一个周期,再做DFT。X(k)只是x(n)的FT在某种程度上的近似, X(k)是x(n)频谱(DTFT)的抽样值。第三章 用DFT对信号进行频谱分析 3.1频谱分析的概述信号的频谱分析就是计算信号的傅里叶变换。连续信号与系统的傅里叶

8、分析显然不能直接用计算机进行计算，使其应用受到限制。而DFT是一种时域和频域均离散化的变换，适合数值计算，成为分析离散信号和系统的有力工具。对连续信号和系统，可以通过时域采样，应用DFT进行频谱分析。用DFT进行频谱分析的基本原理和方法：例如已经知道单位圆上的z变换就是序列的傅里叶变换，即： X（jw）=X（z）|z=jw (3-1)如果对序列x（n）进行N点DFT，得到X（k），X（k）是在区间0,2上对X（jw）的N点等间隔采样。因DFT有FFT算法，故常用DFT对有限长序列进行谱分析。MATLAB程序：syms t;x=sin(t)+cos(t);ezplot(x,-30 30);gri

9、d on;xlabel(' 时间');ylabel('幅值');title(' 时域信号波形');n=0:20; T=0.02;y=sin(n*T)+cos(n*T); figure,subplot(2,2,1);stem(n,y);grid on;xlabel(' n');ylabel('x(n) ');title('500Hz采样信号');Y=fft(y,16);n1=0:15;omega=2*pi/16*(n1-16/2);subplot(2,2,3);stem(omega,abs(fftsh

10、ift(Y);grid on;xlabel('频率');ylabel('幅值');title('n=20时采样信号的频谱');n=0:100; T=0.02;y=sin(n*T)+cos(n*T); subplot(2,2,2);stem(n,y);grid on;xlabel(' n');ylabel(' x(n)');title('500Hz采样信号');Y=fft(y,128);n1=0:127;omega=2*pi/128*(n1-128/2);subplot(2,2,4);stem(ome

11、ga,abs(fftshift(Y);grid on;xlabel('频率');ylabel('幅值');title('n=100时采样信号的频谱');用程序运行后得： 图3-1 时域信号波形 图3-2 500Hz采样点信号和不同采样点信号频谱第四章 吉布斯效应 4.1 吉布斯效应的定义 将具有不连续点的周期函数（如矩形脉冲)进行傅立叶级数展开后，选取有限项进行合成。当选取的项数越多，在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时，该峰起值趋于一个常数，大约等于总跳变值的9%。这种现象称为吉布斯效应。 4.2 吉布斯效应的实

12、现吉布斯现象是当用信号的谐波分量的和来表述具有间断点的波形时出现，并能观察。（1）信号中频率较低的谐波分量的幅值较大，占主体地位，信号波形中所含的频率吉布斯现象越突出。（2）当截取窗变长时，跳变峰向间断点靠近，但跳变峰值并未明显减小，跳变峰所包围的面积减小，通过MATLAB使这种吉布斯现象得到清楚的表现。下面对方波、锯齿波信号进行傅里叶级数展开成三角函数表达式，并利用MATLAB软件与其精准图像进行对比，并通过数值分析验证吉布斯现象。4.2.1方波 周期方波信号: 由于，任何周期函数在满足狄利克雷条件下，可将方波x(t)展成傅里叶级数形式， （4-1）Matlab程序如下：k=inp

13、ut('取点间隔k=');N=input('选取项数N=');t=-2*pi:k:2*pi;s=square(t);y=0;for n=1:2:N, y=y+4/pi*sin(n*t)/n; %前n项和endplot(t,s,t,y);w=(max(y)-1)/2 %峰起值占总跳变值的比例 Matlab分析图像如下所示： k=0.01 N=20 w=8.99% 图4-1 周期方波吉布斯现象k=0.01 N=200 W=8.77% 图4-2周期方波吉布斯现象 k=0.01 N=2000 w=8.94% 图4-3 周期方波吉布斯现象取k(取点间隔)=0.001,0.

14、0001···的情况下，不同N值（选取项数）=200,2000···对应的图像与以上图像类似，故略4.4.2锯齿波周期锯齿波信号：其傅里叶展开表达式为 （4-2）Matlab程序如下：k=input('取点间隔k=');N=input('选取项数N=');t=-2*pi:k:2*pi;s=sawtooth(t-pi);y=0;for n=1:1:N, y=y+2/pi*sin(n*t+(n+1)*pi)/n; %前n项和endplot(t,s,t,y);w=(max(y)-1)/2 %峰起值占总跳变值的

15、比例Matlab分析图像如下所示：k=0.01, y=20, w=6.54% 图4-4 周期锯齿波吉布斯现象 k=0.01, y=200, w=8.55% 图4-5 周期锯齿波吉布斯现象k=0.01 N=2000 w=8.91%图4-6 周期锯齿波吉布斯现象取k(取点间隔)=0.001,0.0001···的情况下，不同N值（选取项数）=200,2000···对应的图像与以上图像类似，故略所以，取点间隔k=0.001逼近度N=200 w =0.0870图4-7 周期锯齿波吉布斯现象第五章 栅栏效应 5.1 栅栏效应的定义快速傅立叶变换得

16、到的频谱是离散谱，是信号的频谱与一个窗函数的频谱做卷积后，按归一化频率分辨率等间隔频域采样的结果，它只给出频谱在离散点上的值，而无法反映这些点之间的频谱内容，即使在其它点上有重要的峰值也会被忽略。这就好像在百页窗内观察窗外的景色，看到的是百叶窗窗缝内的部分景色，而无法看到被百叶窗遮挡住的部分，这就是栅栏现象。5.2 栅栏效应的验证信号包含三种频率成分，分别是f1=20HZ、f2=20.5HZ、f3=40HZ，采样频率fs=100HZ最小记录点数易知是fs/（20.5-20）=200，因此当频域采样点数N>=200时，不出现栅栏效应，而当N<200时，会有栅栏效应误差出现。为了更好的

17、分析DFT栅栏效应，DFT分三种情况：在128点有效数据不补零情况下的分辨率，在128点有效数据且补零至512点情况下分辨率，在512点有效数据下分辨率。5.2.1信号采样频谱分析时给出连续时间信号，需要根据其最高频率确定采样频率，否则将出现频率混叠现象。本设计中在基本谱分析时，采用100Hz的采样率，满足采样定理的要求。故避免了混叠。由公式得采样序列为：x（n）=sin（2f1 n/fs）+sin（2f2n/fs）+sin（2f3/fs）。5.2.2序列FFT设计时虽然给出了参考教材矩阵定义DFT写出的DFT运算函数，但考虑到运算效率，没有调用。而是使用DFT的快速算法FFT。快速傅里叶变换

18、FFT并不是与DFT不相同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小点数DFT的组合，从而减少运算量。常用的FFT是以2为基数，其长度：。它的运算率高，程序比较简单，使用也十分的方便，当需要进行变换的序列长度不是2的整数次方的时候，为了使用以2为基数的FFT,可以用末尾补0的方法，使其长度延伸到2的整数次方。这对实验结果没有影响。又傅里叶变换后的幅值-频率图是以N/2为对称轴的因为如果令A(k)=|X(k)|，可以看出函数自变量是k，也就是说横坐标为0：N-1。但为了让横坐标有物理意义，一般的将其频率0：Fs-Fs/N线性映射到0：N-1

19、（映射关系为f=n*（Fs-Fs/N）/（N-1）即令L(f)=k,A(k)=f(L(f)=|X(L(f)|,如果再令B(f)=|X(L(f)|,即f变成了新函数的自变量（横坐标为频率）。所以Fs/2对应于N/2(即当n=N/2时，f=N/2)，就出现了FFT变换后的频谱图是关于奈奎斯特频率Fs/2对称的。因此我们只选64点和256点对序列进行FFT。5.2.3程序设计 N1=128;N2=512; fs=100;f1=20;f2=20.5;f3=40; n1=0:N1-1;n2=0:N2-1;xn1=sin(2*pi*f1*n1/fs)+sin(2*pi*f2*n1/fs)+sin(2*pi

20、*f3*n1/fs);【在128点有效数据不补零情况下的分辨率】 xk11=fft(xn1,N1)mxk11=abs(xk11(1:N1/2); figure(1); subplot(211);plot(n1,xn1);xlabel('n');title('x(n) 0<=n<127');axis(0,128,-3,3); k1=(0:N1/2-1)*fs/N1; subplot(2，1，2)；plot(k1,mxk11); xlabel('频率 单位Hz');title('X1(k)的幅度谱');【在128点有效数据

21、且补零至512点情况下分辨率】 xn2=xn1,zeros(1,N2-N1); xk12=fft(xn2,N2); mxk12=abs(xk12(1:N2/2); figure(2);subplot(211);plot(n2,xn2);xlabel('n');title('x(n) 0<=n<=511');axis(0,512,-3,3);k2=(0:N2/2-1)*fs/N2;subplot(2，1，2); plot(k2,mxk12);xlabel('频率 单位Hz');title('x1(k)补零后的幅度谱');

22、【在512点有效数据下分辨率】xn3=sin(2*pi*f1*n2/fs)+sin(2*pi*f2*n2/fs)+sin(2*pi*f3*n2/fs); xk2=fft(xn3,N2);mxk3=abs(xk2(1:N2/2); figure(3); subplot(2，1，1);plot(n2,xn3); xlabel('n');title('x(n) 0<=n=511');axis(0,512,-3,3);k3=(0:N2/2-1)*fs/N2; subplot(2，1，2); plot(k3,mxk3);xlabel('频率 单位Hz'

23、;);title('512点有效数据的幅度谱');5.2.4仿真结果与分析Figure1 图5-1 128点有效数据不补零分辨率图由图5-1可看出，在128点有效数据不补零进行FFT情况下，只有40HZ和20.5HZ有幅度值，20HZ的幅度值没有出现，出现了栅栏效应。Figure2 图5-2 128点有效数据且补零至512点分辨率图 由图5-2看出，在128点有效数据且补零至512点进行FFT情况下，相对图5-1得到了高密度的频谱采样，但是并不能得到频谱的更多细节的信息，20HZ的幅度值仍然没有显示，栅栏效应依旧存在。Figure3 图5-3 512有效数据分辨率图由图5-3可看出，在512点有效数据下进行FFT情况下，20HZ，20.5HZ和40HZ的幅度值都能看到。栅栏效应消除。5.2.5依据Matlab仿真结果，可以得出以下结论：（1）利用尾部补零的方法，可得到对采样更密的采样值，即得到高密度的频谱采样，但是并不能得到频谱的