连续时间信号傅里叶级数分析MATLAB课程设计.doc

Matlab应用实践课程设计

绪论
本次课程是通过MATLAB软件来实现数字信号系统里的相关图像和相关仿真的软件。近年来，MATLAB以其强大的矩阵计算和图像视化功能逐渐为国人所知。MATLAB是mathworks公司的软件产品，MATLAB已经成为一个系列产品：MATLAB主包各种工具（toolbox）。功能丰富的工具箱大致分为两类：功能型工具箱和领域型工具箱。功能型工具箱主要用来扩充MATLAB的符号计算功能﹑图形建模仿真功能﹑文字处理功能以及与硬件实时交互功能，能用于多种学科。而领域型工具箱是专业性很强的，如控制工具（control toolbox）﹑信号处理工具箱（signal processing toolbox）等。MATLAB (MATrix LABoratory)具有用法简易、可灵活运用、程式结构强又兼具延展性[3]。
以下为其几个特色： 功能强的数值运算 - 在MATLAB环境中，有超过500种数学、统计、科学及工程方面的函数可使用，函数的标示自然，使得问 题和解答像数学式子一般简单明了，让使用者可全力发挥在解题方面，而非浪费在电脑操作上。
先进的资料视觉化功能 - MATLAB的物件导向图形架构让使用者可执行视觉数据分，并制作高品质的图形，完成科学性或工程 性图文并茂的文章。
高阶但简单的程式环境 - 做为一种直译式的程式语言，MATLAB容许使用者在短时间内写完程式，所花的时间约为用 FORTRAN 或 C 的几分之一，而且不需要编译(compile)及联结 (link) 即能执行，同时包含了更多及更容易使用的内建 功能。
开放及可延伸的架构 - MATLAB容许使用者接触它大多数的数学原使码，检视运算法，更改现存函数，甚至加入自己的函数使 MATLAB成为使用者所须要的环境。
丰富的程式工具箱 - MATLAB的程式工具箱融合了套装前软体的优点，与一个灵活的开放但容易操作之环境，这些工具箱提 供了使用者在特别应用领域所需之许多函数。现有工具箱有：符号运算（利用Maple V的计算核心执行 ）、影像处理、统计分析、讯号处理、神经网路、模拟分析、控制系统、即时控制、系统确认、强建控 制、弧线分析、最佳化、模糊逻辑、mu分析及合成、化学计量分析[4]。

1 MATLAB简介
1.1 MATLAB语言功能
MATLAB功能丰富，可扩展性强。MATLAB软件包括基本部分和专业扩展两大部分的功能。基本部分包括：矩阵的运算和各种变换；代数和超越方程的求解；数据处理和傅立叶变换；数值部分等等，可以充分满足大学理工科本科的计算需要。扩展部分称为工具箱。它实际上是用MATLAB的基本语句辩称的各种子程序集，用于解决某一方面的专门问题，或实现某一类的新算法。
MATLAB 具有以下基本功能
（1）数值计算功能；
（2）符号计算功能；
（3）图形处理及可视化功能；
（3）可视化建模及动态仿真功能[6]。
1.2 MATLAB语言特点
MATLAB 给用户带来的是最直观，最简洁的程序开发环境。它具有以下特点：
（1）语言简洁紧凑，使用方便灵活，库函数极其丰富。MATLAB 程序书写形式自由，利用起丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。
（2）运算符丰富。由于MATLAB 是用C 语言编写的，MATLAB 提供了和C语言几乎一样多的运算符，灵活使用MATLAB 的运算符将使程序变得极为简短。（3）MATLAB 既具有结构化的控制语句（如for 循环，while 循环，break 语句和if 语句），又有面向对象编程的特性。
（4）程序限制不严格，程序设计自由度大。例如，在MATLAB 里，用户无需对矩阵预定义就可使用。
（5）程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。
（6）MATLAB 的图形功能强大。在FORTRAN 和C 语言里，绘图都很不容易，但在MATLAB 里，数据的可视化非常简单。MATLAB 还具有较强的编辑图形界面的能力。
（7）功能强大的工具箱是MATLAB 的另一特色。MATLAB 包含两个部分：
核心部分和各种可选的工具箱。核心部分中有数百个核心内部函数。其工具箱又
分为两类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计
算功能，图示建模仿真功能，文字处理功能以及与硬件实时交互功能，而学科性工具箱是专业性比较强的，如control, toolbox, signl ，proceessing ，toolbox,commumnication toolbox 等[4]。

2 傅里叶级数基本原理概要
2.1 周期信号的傅里叶分解
设有连续时间周期信号，它的周期为T，角频率，且满足狄里赫利条件，则该周期信号可以展开成傅里叶级数，即可表示为一系列不同频率的正弦或复指数信号之和。傅里叶级数有三角形式和指数形式两种。
1. 三角形式的傅里叶级数[2]：
式中系数，称为傅里叶系数，可由下式求得：
[
2. 指数形式的傅里叶级数[2]：

式中系数称为傅里叶复系数，可由下式求得：

周期信号频谱具有三个特点[1]：
（1） 离散性，即谱线是离散的；
（2） 谐波性，即谱线只出现在基波频率的整数倍上；
（3） 收敛性，即谐波的幅度随谐波次数的增高而减小。
周期信号的傅里叶分解用Matlab进行计算时，本质上是对信号进行数值积分运算。在Matlab中有多种进行数值积分运算的方法，我们采用quadl函数，它有两种其调用形式。 (1) y＝quadl(‘func’, a, b)。 其中func是一个字符串，表示被积函数的.m文件名（函数名）；a、b分别表示定积分的下限和上限。
(2) y＝quadl(@myfun, a, b)。其中“@”符号表示取函数的句柄，myfun表示所定义函数的文件名。[6]
2.3 三角形式和指数形式傅里叶级数及各系数间的关系
傅里叶级数的指数形式和三角形式是等价的，其系数可互相转换。表2-1综合了三角
形式和指数形式傅里叶级数及其系数，以及各系数间的关系。
表2-1 周期函数展开为傅里叶级数
形式
指数形式
三角函数形式

展开式

傅里叶系数

系数间的关系

2.2周期信号的频谱
周期信号经过傅里叶分解可表示为一系列正弦或复指数信号之和。为了直观地表示出信号所含各分量的振幅，以频率（或角频率）为横坐标，以各谐波的振幅或虚指数函数的幅度为纵坐标，可画出幅度－频率关系图，称为幅度频谱或幅度谱。类似地，可画出各谐波初相角与频率的关系图，称为相位频谱或相位谱[2]。

3 用MATLAB实现周期信号的傅立叶级数分解与综合
下面以矩形信号为例介绍用MATLAB来实现周期信号的傅立叶技术的分解与综合。
3.1 合成波形与原波形之间的关系
本文使用的连续周期矩形信号是一个周期为4，占空比为50%，幅值为1的矩形波信号，从上面的分析可以得知，这个矩形波信号可以分解成傅立叶级数也是就无数个不同频率的三角波的叠加，用MATLAB软件可以很容易的画出不同个数谐波叠加形成的合成波的形状（源代码见附录），以便将它们与原信号做比较，给理论分析以正确的实验基础。
下图为分别用1到10次谐波叠加形成的合成波与原信号的比较，可以很快地看出，当叠加的谐波数越多是，与原信号的差别就越小，这样就直接证明了傅立叶级数理论的正确性。

图3-1 不同次数的谐波的合成波与原信号的比较

3.2 吉布斯现象
分析傅立叶级数的公式

满足狄里赫利条件的周期函数表示成的傅立叶级数都收敛。狄里赫利条件如下：
1. 在任何周期内，x(t)必须绝对可积；
2. 在任一有限区间中，x(t)只能取有限个最大值或最小值；
3. 在任何有限区间上，x(t)只能有有限个第一类间断点。
所谓的吉布斯现象就是：在x(t)的不可导点上，如果我们只取x(t)等式右边的无穷级数中的有限项作和X(t)，那么X(t)在这些点上会有起伏[1]。
具体现象如下图所示，以下分别为谐波次数为N=50，N=100，N=500合成波的情况。

图3.2 谐波次数为N=50合成波

图3.3 谐波次数为N=100合成波

图3.4 谐波次数为N=500合成波

从上面的图像中可以看出，当N=500的时候，合成波与原来的方波拟合得非常好，但是在不可导的点上，即为x=-1.5,x=-0.5,x=0.5,x=1.5这样的点的时候，合成波会有较大的波动，这就是非常明显的吉布斯现象。

4 用MATLAB实现周期信号的单边频谱及双边频谱。
用周期矩形波信号为例其中周期T=5，脉冲宽度 τ=1。在下图中，τ为一个脉冲的宽度，T为脉冲的周期，它们是方波信号的两个非常重要的参数，它们的大小关系将直接影响方波信号的形状和性质。
以下皆以方波为例来介绍方波各种参数对其频谱的影响以及单边频谱和双边频谱的关系，方波示意图如下：

图4.1 周期矩形脉冲信号
4.1 单边，双边（幅度，相位）频谱及其关系
4.1.1单边，双边（幅度，相位）
如前所述，周期信号可以分解成一系列正弦（余弦）信号或虚指数信号之和，即

其中，
或 幅度和相位
为了直观地表示出信号所含各分量的振幅或，随频率的变化情况，通常以角频率为横坐标，以各次谐波的振幅或虚指数函数的幅度为纵坐标，画出如图4.2和4.3所示的各谐波的振幅或与角频率的关系图，称为周期信号的幅度（振幅）频谱，简称幅度谱。图中每条竖线代表该频率分量的幅度，称为谱线。各谱线顶点连线的曲线（如图中原点所示）称为频谱包络线，它反映了各谐波分量幅度随频率变化的情况。图4.2中幅度谱为单边幅度谱（用绘制的频谱）。图4.3中幅度谱为双边幅度谱（用绘制的频谱）。
类似地，也可画出各谐波初相角与角频率的关系图，如图4.1和4.2中各谐波初相角与角频率的关系图，称为相位频谱，简称相位谱。图4.2中相位谱为单边相位谱。图4.3中相位谱为双边相位谱。如果为实数，那么可用的正负来表示为0或也可把幅度谱和相位谱画在一张图上。
由图可见，周期信号的谱线只出现在频率为等原周期信号频率的整数倍的离散频率上，即周期信号的频谱是离散谱。下面的分析中用的方波的参数为：T=10τ。

图4.2 周期信号的幅度谱和相位谱（单边）

图4.3周期T=5，脉宽tao=1 的矩形脉冲及单边频谱图

4.1.2 单边，双边频谱关系
单边频谱和双边频谱的区别就是求值的范围不同，单边频谱求的是频率大于0的情况，而双边频谱求的是所有频率的情况，即包括频率小于0的情况，这个区别在上面的两张图中可以非常明显地看出来。
4.2以单边幅度频谱为例，研究脉冲宽度与频谱的关系
首先令方波首期T=10。改变脉冲宽度，就是在图4-1中T值不变的情况下，改变的τ值的大小，同时τ必须小于T。在MATLAB软件里可以比较方便地改变这个值。xsqual=@(x)1/2.*(x==-1/2)+1.*(x>-1/2&x<1/2)+1/2.*(x==-1/2);这个语句是控制τ值的，现在的参数是1/2，所的到的是T=10τ的情况，下面是比较三种不同τ值的情况，分别改变1/2为：1/1、1/4，可以得到T=5τ和T=20τ的频谱图像。

图4-4 T=5τ的情况
图4-5 T=10τ的情况

图4-6 T=20τ的情况
容易看出，在T不变的情况下，减小τ值，可以使频谱变得更密集，增大τ值则可以使频谱变得稀疏，因此，需要在不同的情况下选择不同的τ值，才能是系统变得更加符合实际需要。由于周期T相同，因而相邻谱线的间隔相同；脉冲宽度窄，其频谱包络线第一个零点的频率愈高，即信号带宽愈宽，频带内所含的分量愈多。可见，信号的频带宽度与脉冲宽度τ成反比[2]。

4.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&x<1/2)+1/2.*(x==-1/2)固定为1/2，然后分别将a,b值变为：a=-4,b=4和a=-6,b=6，来研究方波周期对其频谱的影响。

图4-7 a=-4,b=4即T=8时的情况

图4-8 a=-5,b=5即T=10时的情况

图4-9 a=-6,b=6即T=12时的情况

通过观察以上三个图像中第一个零点的位置，不难看出：当方波的周期越大，频谱就越密集，周期越小，频谱就越稀疏，其实这点也不难理解。因为τ值不变，改变T值就等于改变了T=ατ中比例系数α的大小。由于周期脉冲信号的时域宽度不变，这时频谱包络线的零点所在位置不变，而当周期增长时，相邻谱线的间隔减少，频谱变密。如果周期无限增长（这时就成为非周期信号），那么，相邻谱线的间隔将趋近于零，周期信号的离散频谱就过渡到非周期信号的连续频谱。随着周期的增长，各谐波分量的幅度也相应减少。脉冲周期T 愈长，谱线间隔愈小，频谱越稠密；反之，则越稀疏。

5用MATLAB实现典型周期信号的频谱
5.1 周期方波脉冲频谱的MATLAB实现
周期方波脉冲信号如图5.1所示，其幅度为1，脉冲宽度‘占空比’：duty=1/2，周期T=5。

图5.1 周期方波脉冲
编写CTFSdbfb.m函数文件，源程序文件见附录程序。调用函数CTFSdbfb.m，即可绘出方波脉冲的双边频谱如5.2，5.3，图所示，其中周期T和占空比duty可变，修改程序即可得到单边频谱如图5.4，5.5所示。

图5.2 周期为T=5，占空比duty=50的方波脉冲双边频谱

图5.3 周期为T=10，占空比duty=80的方波脉冲双边频谱

图5.4 周期为T=5，占空比duty=50的方波脉冲单边频谱

图5.5 周期为T=10，占空比duty=80的方波脉冲单边频谱

由图可以看出，周期方波信号频谱与周期矩形脉冲信号具有相同的规律，由于方波的周期与脉宽比，因此频谱的第一个过零点内只有两根谱线。

5.2 周期三角波脉冲频谱的MATLAB 实现
周期三角波脉冲如图5.6所示，周期T=5，其幅度为1。

图5.6 周期三角波脉冲

MATLAB:
内置有产生三角波的函数sawtooth(t)，其调用格式为:x= sawtooth(t,width):根据width值的不同产生不同形状的三角波,参数width 是0—1 之间的标量,指定在一个周期之间最大值的位置,width是该位置的横坐标和周期的比值.因而,当width=0.5 时产生标准的对称三角波，当width时（可缺省）产生锯齿波。
编写CTFSsjbshbd.m函数文件，源程序文件见附录程序。调用函数CTFSsjbshbd.m，即可绘出方波脉冲的双边频谱如图5.7，5.8所示，其中周期T和width可变，修改程序即可得到单边频谱如图5.9，5.10所示。

图5.7 周期T=5，width=0.5的三角脉冲双边频谱

图5.8 周期T=10，width=1的三角脉冲双边频谱
图5.5 周期三角波双边幅度频谱

图5.9 周期T=5，width=0.5的三角脉冲单边频谱

图5.10 周期T=10，width=1的三角脉冲单边频谱

6 小结即心得体会
通过这次课程设计让我对MATLAB的强大功能有更深的了解和熟悉。大二上下学期和大三下学期我们做过MATLAB的实验。通过以前的一点点基础，我到图书馆里借了六本与MATLAB有关的书。通过在寝室里看书，我对MATLAB有了更近一步的了解和熟悉。
然后对老师给的题目进行分析。连续时间信号傅里叶级数分析及MATLAB实现。课题要求：深入研究连续时间信号傅里叶级数分析的理论知识，利用MATLAB强大的图形处理功能，符号运算功能以及数值计算功能，实现连续时间周期信号频域分析的仿真波形。
课题内容：一、用MATLAB实现周期信号的傅里叶级数分解与综合。以周期矩形波信号为例，绘出包含不同谐波次数的合成波形，观察合成波形与原矩形波形之间的关系及吉布斯现象。二、用MATLAB实现周期信号的单边频谱及双边频谱。以周期矩形波形信号为例：1、绘出单边（幅度，相位）频谱。2、绘出双边（幅度，相位）频谱。3、比较1和2之间的关系。4、以单边幅度频谱为例，研究脉冲宽度与频谱的关系。5、以单边幅度频谱为例，研究脉冲周期与频谱的关系。三、用MATLAB实现典型周期信号的频谱。1、周期方波脉冲信号的单、双边幅度频谱（T、占空比，二个参数可变）。2、周期三角脉冲信号的单、双边幅度频谱（T、width，二个参数可变）。参考相关的书籍写出相应的公式和相关的程序，并且初步分析其正确性。
接下来我用以前安装的MATLAB7.0，输入相关的程序，并且运行它，检查它的正确性，并且修正。在程序的调试过程中遇到很多的问题。有的程序我根本就不知道它错在哪里，也不知道怎么去改正。只有不断的翻书，不断的修改程序，不断的调试。最终总算勉强的完成了老师给我的任务。
在写相关源程序的时候，我还收索了大量的网站，在网上收索了很多关于MATLAB的资料。在这个过程中我发现网上有很多有用的知识。以后应该多注意，充分合理的利用网络，通过网络来学习东西。
这次课程设计业使我明白了在知识的领域里我还有很多很多的不足，并且再一次的深深的体会到理论和实践之间还有很到的差别。在以后的学习中应该多多的注意实践知识的训练和积累。在以后的学习生活中要不断的开拓自己的动手能力，不断的训练自己的动手能力。这次课程设计让我深深的明白了自己以后该做什么，该怎么去做。
这次课程设计课程设计业让我复习了数字信号系统处理里的相关知识。通过这次试验我复习了连续时间信号傅里叶级数，以及单边频谱，双边频谱，典型周期信号的单边，双边频谱分析。通过这次课程设计我还对mathtype数学公式编辑器有了一定的了解，并且会用它编辑公式。对word也有了进一步的掌握。

致 谢
首先感谢学校的老师为我们精心的安排这次课程设计，让我们有机会进行实践锻炼。对所学的知识得到巩固和加强，同时也让我们对MATLAB有了进一步的了解和锻炼。然后还要感谢我们的梁老师平时对我们的耐心教导。同时也感谢在这次实践活动中帮助过我的同学。

参考文献：
[1]刘泉，阙大顺，郭志强.数字信号处理原理与实现[Z]. 北京：电子工业出版社，2009年
[2]刘泉，江雪梅.信号与系统.北京：高等教育出版社，2006年2月
[3]施阳，李俊. MATLAB语言工具箱--ToolBox实用指南. 西安：西北工业大学出版社，1999年4月
[4]贺兴华，周媛媛，王继阳，周晖.MATLAB7.x图像处理. 北京：人民邮电出版社，2006年11月
[5]梁虹. 信号与系统分析及MATLAB实现[Z]. 北京：电子工业出版社，2002年
[6] 罗建军. MATLAB教程[Z]. 北京：电子工业出版社，2005年

附录：
源程序一：
周期信号的傅里叶级数分解与综合的源程序如下：
T=4;tao=2;w=2*pi/T;
a0=quadl(@singrect,-2,2)/T; % 计算a0
N=10; an=zeros(1,N); bn=zeros(1,N);
for k=1:N
an(k)=quadl(@rectcos, -2, 2, [], [], k, w)*2/T;
% 计算an. quadl中的[], []表示以默认精度进行数值积分
% rectcos函数中的后两个参数;
bn(k)=quadl(@rectsin, -2, 2, [], [], k, w)*2/T; %计算bn;
end;
n=1:1:N;
figure(1);
subplot(1,2,1);plot(n, an, '-o');grid on;
subplot(1,2,2);plot(n, bn, '-o');grid on
t=-6:0.01:6;
x=pulstran(t,-8:4:8,'rectpuls',2); ％生成周期矩形脉冲信号
figure(2);subplot(6,2,1);
plot(t, x);
axis([-8,8,-1,2]); grid on;
% 有限项级数逼近
A0=a0;
AN=sqrt(an.^2+bn.^2);
fiN=-atan(bn./an);
subplot(6,2,2); plot(t,A0/2); grid on; %直流项
wave=a0/2;
for k=1:10
wave=wave+an(k)*cos(k*w*t+fiN(k));
subplot(6,2,k+2);plot(t,wave);grid on;
end
程序中“singrect”、“rectcos”和“rectsin”分别为所预先定义的函数文件。其中singrect.m文件为：
function y=singrect(t);
y=(abs(t)<=1); % 定义单个矩形脉冲函数
rectcos.m文件为：
function y=rectcos(t, n, w);
y=(abs(t)<=1). *cos(n*w*t); % 定义了矩形脉冲与余弦函数的乘积
rectsin.m文件为：
function y=rectsin(tn,w);
y=(abs(t)<=1).*sin(n*w*t); % 定义了矩形脉冲与正弦函数的乘积

源程序二：
吉布斯现象源程序：
t=-2:0.001:2;
N=input('N=');
c0=0.5;
fN=c0*ones(1,length(t));
for n=1:2:N
fN=fN+cos(pi*n*t)*sinc(n/2);
end
figure
plot(t,fN)
axis([-2 2 -0.2 1.2])

源程序三：
双边频谱（幅度，相位）分析源程序：
function [A_sym,B_sym]=CTFS2
syms t n y
if nargin<3;Nf=input('pleas Input 所需展开的最高谐波次数: Nf=');end
T=input('pleas Input 信号的周期=');
if nargin<5;Nn=32;end
y=time_fun_s(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));
for k=1:Nf
A_sym(k+1)=double(vpasubs(As,n,k),Nn));
B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));end
if nargout==0
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=time_fun_e
subplot(3,1,2)
stem(k2,abs(S7));
title('连续时间函数周期矩形脉冲的双边幅度谱')
axis([-80,80,0,0.12])
line([-80,80],[0,0])
line([0,0],[0,0.12])
subplot(3,1,3)
stem(k2,angle(S7));
title('连续时间函数周期矩形脉冲的双边相位谱')
axis([-80,80,-4,4])
line([-80,80],[0,0])
line([0,0],[-4,4])
end

function y=time_fun_s(t)
syms a a1
T=input('pleas Input 信号的周期T=');
M=input('周期与脉冲宽度之比M=');
A=1;tao=T/M;a=tao/2;
y1=sym('Heaviside(ta1)')*A;
y=y1-sym('Heaviside(t-a1)')*A;
y=subs(y,a1,a);
y=simple(y);

function x=time_fun_e
T=5;t=-2*T:0.01:2*T;tao=T/5;
x=rectpuls(t,1);
subplot(3,1,1)
plot(t,x)
hold on
x=rectpuls(t-5,1);
plot(t,x)
hold on
x=rectpuls(t+5,1);
plot(t,x)
title('周期为T=5,脉宽tao=1的矩形脉冲')
axis([-10,10,0,1.2])

源程序四：
单边频谱（幅度，相位）分析源程序
a=-5;
b=5;
n=50;
j=sqrt(-1);
tol=1e-6;
T0=b-a;
xsqual=@(x)1/2.*(x==-1/2)+1.*(x>-1/2&x<1/2)+1/2.*(x==-1/2);
out(1)=1/T0.*quad(xsqual,a,b,tol);
xfun=@(x,k,T)xsqual(x).*exp(-j*2*pi*x*k/T);
for i=1:n out(i+1)=1/T0.*quad(xfun,a,b,tol[],i,T0);
end
out1=out(n+1:-1:2);
out1=[conj(out1),out];
angout=angle(out1);
absout=abs(out1);
n1=[-n:n];
subplot(2,1,1)
stem(n1(n+1:2*n+1),absout(n+1:2*n+1));
subplot(2,1,2)
stem(n1(n+1:2*n+1),angout(n+1:2*n+1));

源程序五：
周期方波脉冲频谱双边频谱分析源程序
function [A_sym,B_sym]=CTFS3
syms t n k y
T=10;
if nargin<4;Nf=input('pleas Input 所需展开的最高谐波次数: ');end
if nargin<5;Nn=32;end
y=time_fun_s(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));
for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn));
B_sym(k1)=double(vpa(subs(Bs,n,k),Nn));end
if nargout==0
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)];
subplot(2,1,1)
x=squ_timefun
T=10;t=-2*T:0.01:2*T;
plot(t,x)
title('T=10，占空比为80%的周期方波脉冲')
axis([-20,20,-1,1.2])
line([-20,20],[0,0])
subplot(2,1,2)
stem(k2,abs(S7));
title('连续时间函数周期方波脉冲的双边幅度谱')
axis([-40,40,0,0.6])
end

function y=time_fun_s(t)
syms a a1
T=10;a=T/2;
y1=sym('Heaviside(t)')*2-sym('Heaviside(t-a1)');
y=y1-sym('Heaviside(t+a1)');
y=subs(y,a1,a);
y=simple(y);

function x=squ_timefun
T=10;t=-2*T:0.01:2*T;duty=80;
x=square(t,duty);

源程序六：
周期方波脉冲单边频谱分析源程序
function [A_sym,B_sym]=CTFS31
syms t n k y
T=10;
if nargin<4;Nf=input('pleas Input 所需展开的最高谐波次数: ');end
if nargin<5;Nn=32;end
y=time_fun_s(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));
for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn));
B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));end
if nargout==0
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;
N=Nf*2*pi/T;
k2=0:2*pi/T:N;
subplot(2,1,1)
x=squ_timefun
T=10;t=-2*T:0.01:2*T;
plot(t,x)
title('T=10，占空比为80%的周期方波脉冲')
axis([-20,20,-1,1.2])
line([-20,20],[0,0])
subplot(2,1,2)
stem(k2,abs(S5));
title('周期方波脉冲的单边幅度谱')
axis([0,40,0,0.6])
end

源程序七：
周期三角脉冲双边频谱分析源程序
function [A_sym,B_sym]=CTFS4
syms t n k y
T=10;
if nargin<4;Nf=input('pleas Input 所需展开的最高谐波次数: ');end
if nargin<5;Nn=32;end
y=time_fun_s(t);
A0=2*int(y,t,0,T)/T;
As=int(2*y*cos(2pi*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));
for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k),Nn));
B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn));end
if nargout==0
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)];
subplot(2,1,1)
x=sjb_timefun