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

1、佳木斯大学课程设计报告基于MATLAB的信号吉布斯效应和栅栏效应的验证学 院 信息电子技术专 业 班 级 学 籍 号 姓 名 指导教师 佳 木 斯 大 学2015年7月5日第一章 Matlab和GUI简介1.1 Matlab简介MATLAB（矩阵实验室）是MATrix LABoratory的缩写，是一款由美国The MathWorks公司出品的商业数学软件。MATLAB是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。除了矩阵运算、绘制函数/数据图像等常用功能外，MATLAB还可以用来创建用户界面及与调用其它语言（包括C，C+和FORTRAN）编写的程序。尽管

2、MATLAB主要用于数值运算，但利用为数众多的附加工具箱（Toolbox）它也适合不同领域的应用，例如控制系统设计与分析、图像处理、信号处理与通讯、金融建模和分析等。另外还有一个配套软件包Simulink，提供了一个可视化开发环境，常用于系统模拟、动态/嵌入式系统开发等方面。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本

3、数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C+，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。1.2 GUI简介GUI是Graphical User Interface 图形用户界面的意思，象很多高级编程语言一样，Matlab也

4、有图形用户界面开发环境，随着计算机技术的飞速发展，人与计算机的通信方式也发生的很大的变化，从原来的命令行通讯方式（例如很早的DOS系统）变化到了现在的图形界面下的交互方式，而现在绝大多数的应用程序都是在图形化用户界面下运行的。GUI的广泛应用是当今计算机发展的重大成就之一，他极大地方便了非专业用户的使用。人们从此不再需要死记硬背大量的命令，取而代之的是可以通过窗口、菜单、按键等方式来方便地进行操作。而嵌入式GUI具有下面几个方面的基本要求：轻型、占用资源少、高性能、高可靠性、便于移植、可配置等特点。第二章 离散傅立叶变换2.1 傅立叶变换的定义所谓傅里叶变换就是在以时间为自变量的“信号”与以频

5、率为自变量的“频谱”函数之间的某种变换关系。随时间自变量形式的不同，其傅里叶变换的形式也有不同，常用的两种傅里叶变换：周期序列的离散傅里叶级数(DFS)和非周期序列的傅里叶变换(DTFT)，其表示式分别为： (2-1) (2-2)在实际工作中，当用数字计算机对信号进行频谱分析时，要求信号必须以有限长度的离散值作为输入，而计算所得的频谱值自然也是有限、离散的。上述两种形式的傅里叶变换中，DFS变换满足时、频域自变量的离散化，但其时间变量和频率变量又同时具有周期性；DTFT变换满足时间自变量的有限长度(非周期能量有限信号)，但其频率变量为连续形式。可见，这两种变换都难以实际应用。考虑到DFS变换的

6、时、频域形式虽是周期序列，但每个周期却只有N个独立的复值，知道其一个周期的内容即可得到其它的内容。因此，若从DFS变换的时、频域各取出一个周期，即可构造出时间和频率自变量皆为离散、有限长度的傅氏变换，这就是离散傅里叶变换(DFT)的引出思想，下面进行具体推导。设是一个长度为M的有限长序列，由周期序列与有限长序列的本质联系，可以N()为周期将展开为无重叠的周期序列，即周期延拓为 (2-3)再对进行DFS变换，得到周期离散的频谱，取的主值序列，代入DFS反变换公式即 (2-4)分析可见：在DFS正变换中，只要把一个周期内的乘以对应的，即可得任意k下的；同理，在DFS反变换中，仅用的一个周期的值，即

7、可得到任意n下的。如果同时限制(2-3)式中的n和(2-4)式中的k，使其都只在区间内取值，就得到了一个周期的和一个周期的间的对应关系 (2-5) (2-6)式中，N为DFT变换区间的长度，上两式即称为有限长序列的离散傅里叶变换对。(2-5)式称为离散傅里叶变换，简称DFT；(2-6)式称为离散傅里叶逆变换(Inverse Discrete Fourier Transform)，简称IDFT。由上述DFT的推导过程和定义可以看出，DFT变换具有隐含周期性。一方面由于DFT变换由DFS变换引出，另一方面由DFT定义式中的满足下述周期性 m为整数 (2-7)使得有限长序列和均具有周期性，即满足因此

8、，周期序列可以看作长度为N的有限长序列的周期延拓序列，如式 (2-3)，又可表示为 (2-8)式中是矩形序列，是余数运算表达式，即若 ()，则，表示以N为周期的周期序列。而可看作是从到的主值区间，表示 (2-9)x(n)nn.图2.1有限长序列及其周期延拓上述关系如图2.1所示。同样可将表示为 （2-10） （2-11）上述分析同样说明：有限长序列的离散傅里叶变换，正好是的周期延拓序列的离散傅里叶级数系数的主值序列。也可以理解为，离散傅里叶变换(DFT)是对序列的傅里叶变换(DTFT)的频域采样，运用频域采样方法同样可以推导出DFT的变换公式。傅立叶分析: 建立以时间为自变量的信号和以频率为自

9、变量的频谱函数之间的某种关系,在1822年, 由法国科学家 Fourier提出,基本思想: 任意函数可分解为无穷多个不同频率正弦信号的和, 即频谱分析。 2.2 离散FT(DFT)的用法对常遇到的非周期序列, 有限长或无限, 只能作DTFT, 即连续频谱, X(), 模拟在计算机上做数值计算,实际中, 把N点序列视为一周期序列的一个周期,再做DFT。X(k)只是x(n)的FT在某种程度上的近似, X(k)是x(n)频谱(DTFT)的抽样值。周期信号可以分解成一系列正弦（余弦）信号或虚指数信号之和，即 其中， 或 幅度和相位 为了直观地表示出信号所含各分量的振幅或，随频率的变化情况，通常以角频率

10、为横坐标，以各次谐波的振幅或虚指数函数的幅度为纵坐标，画出如图2和图4所示的各谐波的振幅或与角频率的关系图，称为周期信号的幅度（振幅）频谱，简称幅度谱。图中每条竖线代表该频率分量的幅度，称为谱线。各谱线顶点连线的曲线（如图中原点所示）称为频谱包络线，它反映了各谐波分量幅度随频率变化的情况。类似地，也可画出各谐波初相角与角频率的关系图，如图3和5中各谐波初相角与角频率的关系图，称为相位频谱，简称相位谱。如果为实数，那么可用的正负来表示为0或也可把幅度谱和相位谱画在一张图上。第三章 用DFT对信号进行频谱分析 3.1频谱分析的概述所谓信号的频谱分析就是计算信号的傅里叶变换。连续信号与系统的傅里叶分

11、析显然不能直接用计算机进行计算，使其应用受到限制。而DFT是一种时域和频域均离散化的变换，适合数值计算，成为分析离散信号和系统的有力工具。对连续信号和系统，可以通过时域采样，应用DFT进行频谱分析。用DFT进行频谱分析的基本原理和方法：已经知道单位圆上的z变换就是序列的傅里叶变换，即： X（jw）=X（z）|z=jw如果对序列x（n）进行N点DFT，得到X（k），X（k）是在区间0,2上对X（jw）的N点等间隔采样。因DFT有FFT算法，故常用DFT对有限长序列进行谱分析，实施方法如下：首先依据频率分辨率的要求确定DFT变换区间长度 N 。频谱分析的衡量指标之一是频率分辨率，它是频谱分析中能够

12、分辨的两个相邻频率点谱线的最小间距。在数字频率域，N 点DFT能够实现的频率分辨率是2/N 弧度，进行频谱分析时，要求N>=2/D （D为要求的分辨率）。为了便于使用FFT，一般取 N= 2M. 接下来计算 N 点 DFT ，并以自变量 k 所对应的数字频率 wk=2/N 为为横坐标变量绘制频谱图。用程序运行后得下图：3.2 MATLAB程序：syms t;x=sin(t)+cos(t);ezplot(x,-30 30);grid on;xlabel(' 时间');ylabel('幅值');title(' 时域信号波形');n=0:20;

13、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(fftshift(Y);grid on;xlabel('频率');ylabel('幅值');title('n=20时采样信号的频谱&

14、#39;);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(omega,abs(fftshift(Y);grid on;xlabel('频率');ylabel('幅值');title('n

15、=100时采样信号的频谱');第四章 吉布斯效应 4.1 吉布斯效应的定义 在x(t)的不可导点上，如果我们只取x(t)等式右边的无穷级数中的有限项作和X(t)，那么X(t)在这些点上会有起伏。当选取的项数越多，在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时，该峰起值趋于一个常数，大约等于总跳变值的9%。这种想像成为吉布斯现象。 4.2 吉布斯效应的实现吉布斯现象是当用信号的谐波分量的和来表述具有间断点的波形时出现，并能观察。（1）信号中频率较低的谐波分量的幅值较大，占主体地位，信号波形中所含的频率吉布斯现象越突出。（2）当截取窗变长时，跳变峰向间断点靠近，但跳变

16、峰值并未明显减小，跳变峰所包围的面积减小，通过matlab使这种吉布斯现象得到清楚的表现。模拟信号x（t）=sin(10t)/t 的时域图及其频谱图界面如下：4.3 MATLAB程序：显示波形1a=str2double(get(handles.edit1,'string');c=str2double(get(handles.edit3,'string');syms t;x=sin(c*t)./t;axes(handles.axes1);ezplot(x,-a,a);xlabel('时间');ylabel('幅值');axis(-a

17、 a -4 12);grid on;%hold on;X=fourier(x);XX=abs(X);axes(handles.axes2);ezplot(XX,-3*c,3*c);xlabel('频率');ylabel('幅值');grid on;%hold on;显示波形2b=str2double(get(handles.edit2,'string');c=str2double(get(handles.edit3,'string');ts=-b+1;te=b-1;n=100;t1=linspace(ts,te,n);x1=sin

18、(c*t1)./t1;axes(handles.axes3);plot(t1,x1)axis(-b,b,-4,12);xlabel('时间');ylabel('幅值');grid on;%hold on;X1=abs(fft(x1,128);fs=n/(te-ts);f=(0:length(X1)/2-1)*fs/128*2*pi;axes(handles.axes4);plot(f,X1(1:length(X1)/2);xlabel('频率');ylabel('幅值');grid on;%hold on;第五章 栅栏效应 5.1

19、 栅栏效应的定义快速傅立叶变换得到的频谱是离散谱，是信号的频谱与一个窗函数的频谱做卷积后，按归一化频率分辨率等间隔频域采样的结果，它只给出频谱在离散点上的值，而无法反映这些点之间的频谱内容，即使在其它点上有重要的峰值也会被忽略。这就好像在百页窗内观察窗外的景色，看到的是百叶窗窗缝内的部分景色，而无法看到被百叶窗遮挡住的部分，这就是栅栏现象。5.2 栅栏效应的验证设计一个离散傅立叶计算程序，计算方波的频谱。观察栅栏效应带来的计算误差。设计如下：利用505Hz方波信号的频谱分析说明栅栏效应所造成的频谱误差。设定采样频率：fs=5120，MATLAB中默认的FFT计算点数为512，其离散频率点为：f

20、i=i*fs/N=i*5120/512=10,i=0,1,2,3,N/2位于505HZ位置的真实谱峰被挡住看不见，看见的只是它们在相邻频率500Hz或510Hz处能量泄漏的值。若设fs=2560Hz,则频率间隔为df=5Hz，重复上述分析步骤，这时在505 Hz位置有谱线，因此可以得到它们的精确值。从时域看，这个条件相当于对信号进行整周期采样，实际中常用此方法来提高周期信号的频谱分析精度。将上图中的后两个图是不同采样频率下的采样图，从图看出是不一样的，下面我们把后面的两图，即刚才离散的采样的频率图进行局部放大可以得到下面的图，由此看出栅栏效应的影响。我们知道，N点DFT是频率区间上对信号的频谱

21、进行N点等间隔采样，而采样点之间的频谱函数值是不知道的。这就好像从(N+1)个栅栏缝隙中观看信号的频谱情况，仅得到N个缝隙中看到的频谱函数值。因此称这种现象为栅栏效应。由于栅栏效应，有可能漏掉(挡住)大的频谱分量。为了把原来被“栅栏”挡住的频谱分量检测出来，需要提高频率分辨率(用频率采样间隔F来描述，表示谱分析中能够分辨的两个频率分量的最小间隔)，也就是说：间隔F越小，谱分析越接近原连续频谱，频率分辨率越高。按照(1-10)式，为使F减小，需要减小采样频率或增大采样点数N。由于采样频率的降低会引起谱分析范围的减少或者引起频谱混叠，所以在维持不变的条件下。频率分辨率可以增加采样点数N，因为，只有

22、增加对信号的观察时间，才能增加N。和N可以按照下面两式进行选择： (5-1) (5-2)增加观察时间或者采样点数最简单的办法可以采用在原序列尾部补零的方法，改变序列长度N(即改变DFT变换区间长度)，从而增加了频域采样点数和采样点位置，使原先漏掉的某些频谱分量被检测出来。界面图5.3 matlab程序采样程序如下：t1=0:1/5120:0.1;m1=0:10:5120;t2=0:1/2560:0.2;m2=0:5:2560;y1=square(2*pi*505*t1);transf1=abs(fft(y1)/256);y2=square(2*pi*505*t2);transf2=abs(ff

23、t(y2)/256);subplot(2,2,1);plot(t1,y1);title('square(2*pi*505*t1)');subplot(2,2,3);stem(m1(1:256);transf1(1:256);title('fs=5120');subplot(2,2,2);plot(t2,y2);title('square(2*pi*505*t2)');subplot(2,2,4);stem(m2(1:256),transf2(1:256);title('fs=2560');界面显示程序如下：显示波形1a=str2d

24、ouble(get(handles.edit4,'String');t1=0:1/a:0.1;m1=0:10:a;y1=square(2*pi*505*t1);transf1=abs(fft(y1)/256);axes(handles.axes1);stem(t1,y1);axes(handles.axes2);stem(t1,transf1);axes(handles.axes7);stem(t1,y1);axes(handles.axes8);stem(t1,transf1);axis(0.0525 0.0575 0 0.5);显示波形2b=str2double(get(handles.edit3,'String');t2=0:1/b:0.1;m2=0:5:b;y2=square(2*pi*505*t2);transf2=abs(fft(y2)/256);axes(handles.axes3);stem(t2,y2);axes(handles.axes4);stem(t2,transf2);axes(handles.axes5);stem(t2,y2);axes(handles.axes