连续时间信号傅里叶级数分析报告及matlab实现文档格式.docx
《连续时间信号傅里叶级数分析报告及matlab实现文档格式.docx》由会员分享,可在线阅读,更多相关《连续时间信号傅里叶级数分析报告及matlab实现文档格式.docx(29页珍藏版)》请在冰豆网上搜索。
傅里叶级数;
周期信号;
频谱
绪论
在科学技术飞速发展的今天,计算机正逐步将科技人员从繁重的计算工作中解脱出来。
在进行科学研究与工程应用中,往往需要大量的科学计算,一些科技人员曾经尝试使用传统的高级语言Basic、Fortran及C语言编写程序,以减轻工作量。
但编制程序需要掌握高级语言的语法,还要对各种算法进行了解,这对大多数科技人员来说是不大现实的,而且也是没有没有必要的。
MATLAB正是在这一应用要求背景下产生的数学类科技应用软件。
它具有的顶尖的数值计算功能、强大的图形可视化功能及简洁易学的“科学便捷式”工作环境和编程语言,从根本上满足了科技人员对工程数学计算的要求,并将科技人员从繁重的数学运算中解放出来,因而越来越受到广大科技工作者的普遍欢迎[1]。
MATLAB是matrix和laboratory前三个字母的缩写,意思是“矩阵实验室”,是MathWorks公司推出的数学类科技应用软件。
其Dos版本(MATLAB1.0)发行于1984年,现已推出了Windows版本(MATLAB5.3)。
经过十多年的不断发展与完善,MATLAB已发展成为由MATLAB语言、MATLAB工作环境、MATLAB图形处理系统、MATLAB数学函数库和MATLAB应用程序接口五大部分组成的集数值计算、图形处理、程序开发为一体的功能强大的系统。
MATLAB由“主包”和三十多个扩展功能和应用学科性的工具箱(Toolboxs)组成。
目前,MATLAB已经成为国际上最流行的电子仿真计算机辅助设计的软件工具,现在的MATLAB已经不仅仅是一个“矩阵实验室(MatrixLaboratory)”,它已经成为一种实用的、全新的计算机高级语言。
正是由于MATLAB在数值计算及符号计算等方面的强大功能,使MATLAB一路领先,成为数学类科技应用软件中的佼佼者。
目前,MATLAB已成为国际上公认的最优秀的科技应用软件。
MATLAB的上述特点,使它深受工程技术人员及科技专家的欢迎,并很快成为应用学科计算机辅助分析、设计、仿真、教学等领域不可缺少的基础软件。
1MATLAB简介
1.1MATLAB语言功能
MATLAB是一个高精度的科学计算语言,它将计算、可视化编程结合在一个容易使用的环境中,在这个环境中,用户可以把提出的问题和解决问题的办法用熟悉的数学符号表示出来,它的典型使用包括:
(1)数学和计算;
(2)运算法则;
(3)建模、仿真;
(4)数值分析、研究和可视化;
(5)科学的工程图形;
(6)应用程序开发,包括创建图形用户接口。
1.2MATLAB语言特点
MATLAB是一个交互式系统,他的基本数据单元是数组,这个数组不要求固定的大小,因此可以让用户解决许多技术上的问题,特别是那些包含矩阵和矢量运算的问题。
MATLAB的指令表达与数学、工程中常用的习惯形式相似,与C、Fortran、等高级语言相比,它的语法规则更简单、表达更符合工程习惯,正因为如此,人们用MATLAB语言编写程序就犹如在便笺上书写公式和求解,因而MATLAB被称为“便笺式”的科学工程语言。
MATLAB的最重要特征使他拥有解决特定应用问题的程序组,也就是TOOLBOX(工具箱),如信号处理工具箱,控制系统工具箱、神经网络工具箱、模糊逻辑工具箱、通信工具箱和数据采集工具箱等许多专用工具箱,对大多数用户来说,要想灵活、高效地运用这些工具箱,通常都需要学习相应的专业知识。
此外,开放性也许是MATLA最重要和最受欢迎的特点之一。
除内部函数外,所有的MATLAB主要文件和各工具箱文件都是可读的、可改的源文件,因为工具箱实际上是有一组复杂的MATLAB函数(M文件)组成,它扩展了MATLAB的功能,用以解决待定的问题,因此用户可以通过对源文件进行修改和加入自己编写的文件去构建新的专用工具箱。
2连续时间周期信号的傅里叶级数
频域分析法即傅里叶分析法,它是变换域分析法的基石。
其中,傅里叶级数是变换域分析法的理论基础,傅里叶变换作为频域分析法的重要数学工具,具有明确的物理意义,在不同的领域得到广泛的应用
2.1连续时间周期信号的分解
以高等数学的知识,任何周期为T的周期函数
,在满足狄里赫利条件时,则该周期信号可以展开成傅里叶级数。
傅里叶级数有三角形式和指数形式两种[2]。
2.1.1三角形式的傅里叶级数
三角形式的傅里叶级数为:
式中系数
、
称为傅里叶系数,可由下式求得。
其中,
为基波频率,
为n次谐波频率。
如果将
式中同频率的正弦和余弦分量合并,则三角形式的傅里叶级数可表示为:
上式中
可以看出,傅里叶系数
和
都是
或
的函数,其中
是
的偶函数,即有
;
而
的奇函数,即有
。
2.1.2指数形式的傅里叶级数
根据欧拉公式:
并考虑
奇偶性可将
改写为指数形式的傅里叶级数:
即周期信号可分解为一系列不同频率的虚指数信号之和,式中
称为傅里叶复系数,可由下式求得:
2.2连续时间周期信号的傅里叶综合
任何满足狄里赫里条件的周期信号,可以表示成式
的和式形式,
式常称为连续周期信号的傅里叶级数综合公式。
一般来说,傅里叶级数系数有无限个非零值,即任何具有有限个间断点的周期信号都一定有一个无限项非零系数的傅里叶级数表示。
但对数值计算来说,这是无法实现的。
在实际的应用中,但我们可以用有限项的傅里叶级数求和来逼近。
为了比较有限项谐波的逼近情况,本次课设编写了程序来绘制波形以给读者一个直观的感受。
调用xiebo.m函数文件,即可绘出周期矩形波信号各次谐波的合成波形。
如图2.1所示。
图2.1周期矩形脉冲信号的合成
由图2.1可见,当它所包含的谐波分量越多时,合成波形愈接近于原来的矩形波脉冲(。
由图2.1还可以看到,合成波形所包含的谐波分量愈多时,除间断点附近外,它越接近于原矩形波脉冲。
在间断点附近,随着所含谐波次数的增加,合成波形的尖峰愈接近间断点,但尖峰幅度并未明显减少。
可以证明,即使合成波形所含谐波次数
时,在间断点处仍有约9%的偏差,这种现象称为吉布斯(Gibbs)现象。
在傅里叶级数的项数取得很大时,间断点处尖峰下的面积非常小以致趋近于零,因而在均方的意义上合成波形同原波形的真值之间没有区别[4]。
2.3吉布斯现象
上一节中我们提到了吉布斯现象,本节我们将作重点来讨论。
我们知道满足狄里赫利条件的周期函数表示成的傅立叶级数都收敛。
狄里赫利条件如下:
1.在任何周期内,x(t)必须绝对可积;
2.在任一有限区间中,x(t)只能取有限个最大值或最小值;
3.在任何有限区间上,x(t)只能有有限个第一类间断点。
所谓的吉布斯现象就是:
在x(t)的不可导点上,如果我们只取x(t)等式右边的无穷级数中的有限项作和X(t),那么X(t)在这些点上会有起伏[1]。
具体现象如下图所示,以下分别为谐波次数为N=50,N=100,N=500合成波的情况。
图2.2不同时N值时的合成波
从上面的图像中可以看出,当N=500的时候,合成波与原来的方波拟合得非常好,但是在不可导的点上,即为x=-1.5,x=-0.5,x=0.5,x=1.5这样的点的时候,合成波会有较大的波动,这就是非常明显的吉布斯现象。
3连续时间周期信号的频谱分析
3.1单边与双边频谱关系
如前所述,周期信号可以分解成一系列正弦(余弦)信号或虚指数信号之和,为了直观地表示出信号所含各分量的振幅
,随频率的变化情况,通常以角频率为横坐标,以各次谐波的振幅
或虚指数函数
的幅度为纵坐标,画出如图3.1和3.2所示的各谐波的振幅
与角频率的关系图,称为周期信号的幅度(振幅)频谱,简称幅度谱。
图中每条竖线代表该频率分量的幅度,称为谱线。
各谱线顶点连线的曲线(如图中原点所示)称为频谱包络线,它反映了各谐波分量幅度随频率变化的情况。
图3.1中幅度谱为单边幅度谱(用
绘制的频谱)。
图3.2中幅度谱为双边幅度谱(用
类似地,也可画出各谐波初相角
与角频率的关系图,如图3.1和3.2中各谐波初相角
与角频率的关系图,称为相位频谱,简称相位谱。
图3.1中相位谱为单边相位谱。
图3.2中相位谱为双边相位谱。
如果
为实数,那么可用
的正负来表示
为0或
也可把幅度谱和相位谱画在一张图上。
由图可见,周期信号的谱线只出现在频率为
等原周期信号频率的整数倍的离散频率上,即周期信号的频谱是离散谱。
图3.1周期信号的单边幅度谱和相位谱
图3.2周期信号的双边幅度谱和相位谱
由此可见周期信号频谱具有三个特点:
(1)离散性,即谱线是离散的;
(2)谐波性,即谱线只出现在基波频率的整数倍上;
(3)收敛性,即谐波的幅度随谐波次数的增高而减小[3]。
单边频谱和双边频谱的区别就是求值的范围不同,单边频谱求的是频率大于0的情况,而双边频谱求的是所有频率的情况,即包括频率小于0的情况,这个区别在上面的两张图中可以非常明显地看出来。
3.2以单边幅度频谱为例,研究脉冲宽度与频谱的关系
首先令方波首期T=5。
改变脉冲宽度,就是在图3.3中T值不变的情况下,改变的τ值的大小,同时τ必须小于T。
在MATLAB软件里可以比较方便地改变这个值。
xsqual=@(x)1/2.*(x==-1/2)+1.*(x>
-1/2&
x<
1/2)+1/2.*(x==-1/2);
这个语句是控制τ值的,现在的参数是1/2。
下面是比较三种不同τ值的矩形脉冲单边频谱图像。
图3.3不同τ值时的频谱
容易看出,在T不变的情况下,减小τ值,可以使频谱变得更密集,增大τ值则可以使频谱变得稀疏,因此,需要在不同的情况下选择不同的τ值,才能是系统变得更加符合实际需要。
由于周期T相同,因而相邻谱线的间隔相同;
脉冲宽度窄,其频谱包络线第一个零点的频率愈高,即信号带宽愈宽,频带内所含的分量愈多。
可见,信号的频带宽度与脉冲宽度τ成反比[2]。
3.3以单边幅度频谱为例,研究脉冲周期与频谱的关系
上面是改变τ值来观察频谱的变化情况,现在来改变T值以达到改变频谱的目的。
在MATLAB代码中a=-5;
b=5;
T0=b-a;
这几句代码是用来控制方波的T值的,b-a就是方波的周期T,在上面的讨论中使用的参数是a=-5,b=5,现在将τ的参数,即这句xsqual=@(x)1/2.*(x==-1/2)+1.*(x>
1/2)+1/2.*(x==-1/2)固定为1/2,然后分别将a,b值变为:
a=-4,b=4和a=-6,b=6,来研究方波周期对其频谱的影响。
图3.4不同T值的频谱
通过观察以上三个图像中第一个零点的位置,不难看出:
当方波的周期越大,频谱就越密集,周期越小,频谱就越稀疏,其实这点也不难理解。
因为τ值不变,改变T值就等于改变了T=ατ中比例系数α的大小。
由于周期脉冲信号的时域宽度不变,这时频谱包络线的零点所在位置不变,而当周期增长时,相邻谱线的间隔减少,频谱变密。
如果周期无限增长(这时就成为非周期信号),那么,相邻谱线的间隔将趋近于零,周期信号的离散频谱就过渡到非周期信号的连续频谱。
随着周期的增长,各谐波分量的幅度也相应减少。
脉冲周期T愈大,谱线间隔愈小,频谱越稠密;
反之,则越稀疏。
4典型周期脉冲的频谱
4.1周期方波脉冲频谱的MATLAB实现
周期方波脉冲信号如图4.1所示,其幅度为1,脉冲宽度‘占空比’:
duty=0.5,周期T=5。
图4.1周期方波脉冲
编写fangbo.m函数文件,源程序文件见附录程序四。
调用函数fangbo.m,即可绘出方波脉冲的双边频谱,其中周期T和占空比duty可变,修改程序即可得到单边频谱。
将在下一小节中给出不同参数时的频谱图。
4.1.1周期方波脉冲双边频谱的MATLAB实现
图4.2周期方波脉冲的双边频谱a
图4.3周期方波脉冲的双边频谱b
图4.4周期方波脉冲的双边频谱c
由上面三个图可以看出,当T一定时占空比越大频谱主瓣的宽度越大,当占空比一定时周期越小频谱的主瓣宽度越大。
周期方波信号频谱与周期矩形脉冲信号具有相同的规律,这里不再赘述。
4.1.2周期方波脉冲单边频谱的MATLAB实现
图4.5周期方波脉冲的单边频谱a
图4.6周期方波脉冲的单边频谱b
图4.7周期方波脉冲的单边频谱c
单边频谱就是双边频谱正半轴部分,其具有的规律也与双边频谱相同。
4.2周期三角波脉冲频谱的MATLAB实现
周期三角波脉冲如图4.8所示,周期T=5,其幅度为1。
图4.8周期三角波脉冲
MATLAB内置有产生三角波的函数sawtooth(t),其调用格式为:
x=sawtooth(t,width);
根据width值的不同产生不同形状的三角波,参数width是0—1之间的标量,指定在一个周期之间最大值的位置,width是该位置的横坐标和周期的比值.因而,当width=0.5时产生标准的对称三角波,当width不等于0.5时(可缺省)产生锯齿波。
在附录的源程序五中,只需要给出不同的T值和width值,就会得到三角波的不同的双边和单边频谱图。
4.2.1周期三角波双边频谱的MATLAB实现
图4.9周期三角波脉冲的双边频谱a
图4.10周期三角波脉冲的双边频谱b
图4.11周期三角波脉冲的双边频谱c
由上面三个图可以看出,当三角波为脉宽width=0.5的对称三角波时频谱图在0点的幅值为零。
主瓣宽度与周期和脉冲宽度的关系与方波时的规律基本一致,这里不再赘述。
4.2.2周期三角波单边频谱的MATLAB实现
图4.12周期三角波脉冲的单边频谱a
图4.13周期三角波脉冲的单边频谱b
图4.14周期三角波脉冲的单边频谱c
5小结即心得体会
本次课程设计至此已经接近尾声,一周的时间虽然很短暂,但在这一个星期的设计过程中收获颇多。
设计的核心内容就是利用MATLAB强大的图形处理功能,符号运算功能以及数值计算功能,实现连续时间周期信号频域分析的仿真波形。
整个设计过程中首先对所学的信号与系统与数字信号处理有了更深的了解,比如傅立叶级数、信号频谱等;
其次,实现过程是通过MATLAB软件完成的,MATLAB的图形功能强大,具有良好的人机界面,此次设计过程中熟练了MATLAB的编程,掌握了很多函数的作用及使用方法;
最后,通过此次课程设计,我对设计所用到的软件MATLAB有了更加深刻地了解,MATLAB不管在数值计算方面的功能很强大,而且其图形仿真功能更能满足各个领域的需要,因此我们以后更要经常运用MATLAB软件,使其成为自己不可或缺的工具。
在写相关源程序的时候,我还收索了大量的网站,在网上收索了很多关于MATLAB的资料。
在这个过程中我发现网上有很多有用的知识。
以后应该多注意,充分合理的利用网络,通过网络来学习东西。
在收集资料的阶段我复习了数字信号系统处理里的相关知识。
对以前的理论知识有了更进一步的认识和理解。
通过这次课程设计我还对mathtype数学公式编辑器有了一定的了解,并且会用它编辑公式。
对word也有了进一步的掌握。
虽然我顺利完成了课程设计的要求,但是我感觉到我对MATLAB的理解我掌握还停留在比较浅的层次。
要想真正掌握它还需要继续努力学习它。
这次课程设计也使我明白了在知识的领域里我还有很多很多的不足,并且再一次的深深的体会到理论和实践之间还有很到的差别。
在以后的学习中应该多多的注意实践知识的训练和积累。
在以后的学习生活中要不断的开拓自己的动手能力,不断的训练自己的动手能力。
这次课程设计让我深深的明白了自己以后该做什么,该怎么去做。
致谢
感谢学校给我们这次MATLAB课程设计的机会,不仅让我们更加学会了MATLAB的强大图形处理方法,掌握了MATLAB的编程技术,而且也锻炼了我们的动手能力。
通过这次课设让我明白了理论联系实践的重要性,书本上的理论知识学了不少,我们必须得应用到实践当中,做到学以致用,这样我们才能有不断的创新。
这次课程设计也感谢指导老师在设计过程中的辅导以及同学们的帮助。
没有他们的帮助我不会那么快克服那些困难,也不会这么快学到这么多的知识。
附录
注意:
由于大部分程序都很相似,这里只给出五个主要的源程序。
对于那些只需要该参数就可以实现的程序也给只出主要程序。
源程序一:
(连续周期信号的分解与综合——谐波分析)
function[A_sym,B_sym]=CTFS1
symstnkx
T=5;
tao=0.2*T;
a=0.5;
ifnargin<
4;
Nf=6;
end
5;
Nn=32;
x=time_fun_x(t);
A0=2*int(x,t,-a,T-a)/T;
%求傅里叶级数展开式的系数。
As=int(2*x*cos(2*pi*n*t/T)/T,t,-a,T-a);
Bs=int(2*x*sin(2*pi*n*t/T)/T,t,-a,T-a);
A_sym
(1)=vpa(A0,Nn);
fork=1:
Nf%Nf为最高次谐波的次数
A_sym(k+1)=vpa(subs(As,n,k),Nn);
B_sym(k+1)=vpa(subs(Bs,n,k),Nn);
ifnargout==0
c=A_sym;
disp(c)
d=B_sym;
disp(d)
t=-8*a:
0.01:
T-a;
f1=2*(0.2/2+0.1871.*cos(2*pi*1*t/5)+0.*sin(2*pi*1*t/5));
f2=2*(0.1514.*cos(2*pi*2*t/5)+0.*sin(2*pi*2*t/5));
f3=2*(0.1009.*cos(2*pi*3*t/5)+0.*sin(2*pi*3*t/5));
f4=2*(0.0468.*cos(2*pi*4*t/5)+0.*sin(2*pi*4*t/5));
f5=2*(-0.0312.*cos(2*pi*6*t/5)+0.*sin(2*pi*6*t/5));
f6=f1+f2;
f7=f6+f3;
f8=f7+f4+f5;
subplot(2,2,1)
plot(t,f1),holdon
title('
基波'
)
subplot(2,2,2)
plot(t,f6),holdon
基波+2次谐波'
subplot(2,2,3)
plot(t,f7),holdon
基波+2次谐波+3次谐波'
subplot(2,2,4)
plot(t,f8),holdon
基波+2次谐波+3次谐波+6次谐波'
functiony=time_fun_e%定义矩形脉冲函数(绘图用)
h=1;
t=-8*a:
e1=1/2+1/2.*sign(t+tao/2);
e2=1/2+1/2.*sign(t-tao/2);
y=h.*(e1-e2);
functionx=time_fun_x(t)%定义矩形脉冲函数(积分用)
x1=sym('
Heaviside(t+0.5)'
)*h;
x=x1-sym('
Heaviside(t-0.5)'
源程序二:
(吉布斯效应)
t=-2:
0.001:
2;
N=input('
N='
);
%N为输入要达到的最高次谐波的次数
c0=0.5;
fN=c0*ones(1,length(t));
forn=1:
2:
N
fN=fN+cos(pi*n*t)*sinc(n/2);
figure
plot(t,fN)
axis([-22-0.21.2])
源程序三:
(矩形脉冲的频谱分析)
function[A_sym,B_sym]=CTFS2
symstny
3;
Nf=input('
pleasInput所需展开的最高谐波次数:
Nf='
T=input('
pleasInput信号的周期T='
y=fun_in(t);
A0=2*int(y,t,0,T)/T;
As=int(2*y*cos(2*pi*n*t/T)/T,t,0,T);
Bs=int(2*y*sin(2*pi*n*t/T)/T,t,0,T);
A_sym
(1)=double(vpa(A0,Nn));
Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn));
B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));
S1=fliplr(A_sym)
S1(1,k+1)=A_sym
(1)
S2=fliplr(1/2*S1)
S3=fliplr(1/2*B_sym)
S3(1,k+1)=0
S4=fliplr(S3)
S5=S2-i*S4;
S6=fliplr(S5);
N=Nf*2*pi/T;
k2=-N:
2*pi/T:
N;
S7=[S6,S5(2:
end)];
x=fun_mc
subplot(3,1,2)
st