x1=sym('Heaviside(t+0.5)')*h; x=x1-sym('Heaviside(t-0.5)')*h;
%------------------------------------------- function y=time_fun_e(t)
% 该函数是CTFShchsym.m的子函它由符号函数和表达式写成 a=0.5;T=5;h=1;tao=0.2*T;t=-8*a:0.01:T-a; e1=1/2+1/2.*sign(t+tao/2); e2=1/2+1/2.*sign(t-tao/2);
y=h.*(e1-e2); %连续时间函数-周期矩形脉冲 （d）程序运行结果
在MATLAB命令窗口键入CTFShchsym，并按回车键，即可绘制出周期矩形波形各次谐波合成波形，如图4-1所示。
图4-1 周期矩形脉冲傅立叶分解与综合结果
执行如下命令就可以得到三角级数展开系数： [A_sym,B_sym]=CTFShchsym 结果为 A_sym =
0.2000 0.3742 0.3027 0.2018 0.0935 0.0000 -0.0624
B_sym =
0 0 0 0 0 0 0 4、周期信号频谱分析
对周期为T1 的信号f(t)进行傅立叶级数展开可得到
f(t)??Fnen????jn?1t?a0??ancos(n?1t)??bnsin(n?1t)
n?1n?1??其中
25
??12? T11(an?jbn) （3-8） 2Fn?Fnej?n?如果求出an和bn，根据以下两式可以画出周期信号的幅度谱Fn~?和相位谱?n~?。
Fn?12an?bn2 幅频 2???arctan 相频
nbnan5、MATLAB实现
由上述分析可知，只要求出了周期信号的傅立叶级数幅度Fn和相位?n，就可以根据Fn和?n随频率??n?1的变化关系画出幅度谱和相位谱。
现仍以周期矩形脉冲为f(t)??G?(t?nT)，（??1,T?5）为例，来说明如
n????何用MATLAB绘制周期信号的频谱图，并对周期信号的频谱特性进行分析。 由于绘制频谱图的前提是必须先求出周期信号的傅立叶级数系数，因此只需对上述给出的求周期信号傅立叶级数的函数CTFShchsym.m进行适当修改，即可编写出绘制周期信号频谱的通用函数。
需注意的是，由于周期信号的频谱是离散的，故在绘制频谱时，采用的是stem命令而不是plot命令，下面给出实现上述通用程序CTFStpshsym.m的过程。 （1）实现流程
此处采用三角形式傅立叶级数分解，求出分解系数an和bn，再利用式（3-8）求出傅立叶复指数系数Fn，画出Fn的幅度谱Fn~?。谐波的阶数Nf可任意指定，此处指定Nf?60。
（a）实现流程如下：
? 编写子函数y=time_fun_s(t)，用符号表达式表示出周期信号在第一个周
期内的符号表达式，并赋值返回给符号变量y； ? 编写子函数x=time_fun_e(t)，求出该周期信号在绘图区间内的信号样值，
并赋值给返回变量x；
? 编写求解信号傅立叶复指数系数Fn的通用函数，该函数流程如下：
① 调用函数time_fun_s(t)，获取周期信号的符号表达式； ② 求出信号的三角形式的傅立叶系数an和bn；
26
③ 求出信号的傅立复指数系数Fn； ④ 绘制Fn的幅度频谱图；
⑤ 调用函数time_fun_e(t)绘制信号波形。
（b）算法提示
? 采用符号积分int求[0,T]内时间信号的三角级数展开系数：a0?A0，
an?As，bn?Bs;
? 用循环语句for…end求出三角级数展开系数an和bn的数值，并赋值给变
量A_sym(k+1)， B_sym(k+1)；
? 从三角级数展开系数an和bn得到复指数展开系数Fn：（3-8）Fn可以根据式
求出，需要注意的是an、bn和Fn的自变量取值情况，an、bn的变量n的取值范围为n?0,1,2,3,?,N，而Fn的变量n的取值范围为
n?0,?1,?2,?3,?,?N，为了从an和bn得
到Fn，需要用到MATLAB的反折函数fliplr来实现频谱的反折；例如已知一个
序列a=[5 4 3 2 1]，采样位置为：0，1，2，3，4。要组成员序列反折后再叠加原序列的一个新序列（采样位置：-4，-3，-2，-1，0，1，2，3，4）。程序如下：n=0:4;
a=[5 4 3 2 1];
subplot(2,1,1);stem(n,a); b=fliplr(a); k=-4:4;
c=[b,a(2:end)];
subplot(2,1,2);stem(k,c);
结果如右图所示。
（c）源程序代码
编写函数文件CTFStpshsym.m，如下所示： [CTFStpshsym.m]
function [A_sym,B_sym]=CTFSshbpsym(T,Nf)
% 采用符号计算求[0,T]内时间函数的三角级数展开系数。 % 函数的输入输出都是数值量 % Nn 输出数据的准确位数
27
% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数
% B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % T T=m*tao， 信号周期 % Nf 谐波的阶数 % Nn 输出数据的准确位数
% m (m=T/tao)周期与脉冲宽度之比,如m=4,8,16,100等 % tao 脉宽:tao=T/m syms t n y
if nargin<3;Nf=input('please Input 所需展开的最高谐波次数:Nf=');end T=input('please Input 信号的周期T='); 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); %对A_sym阵左右对称交换
S1(1,k+1)=A_sym(1); %A_sym的1*k阵扩展为1*(k+1)阵
S2=fliplr(1/2*S1); %对扩展后的S1阵左右对称交换回原位置 S3=fliplr(1/2*B_sym); %对B_sym阵左右对称交换 S3(1,k+1)=0; %B_sym的1*k阵扩展为1*(k+1)阵
S4=fliplr(S3); %对扩展后的S3阵左右对称交换回原位置
S5=S2-i*S4; % 用三角函数展开系数A、B值合成付里叶指数系数 S6=fliplr(S5); N=Nf*2*pi/T; k2=-N:2*pi/T:N; S7=[S6,S5(2:end)]; subplot(3,3,3)
x=time_fun_e(t) % 调用连续时间函数-周期矩形脉冲 subplot(2,1,1)
stem(k2,abs(S7)); %画出周期矩形脉冲的频谱（T=M*tao） title('连续时间函数周期矩形脉冲的双边幅度谱') axis([-80,80,0,0.12]) line([-80,80],[0,0]) line([0,0],[0,0.12]) end
28
%------------------------------------------- function y=time_fun_s(t)
% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 syms a a1
T=input('please Input 信号的周期T='); M=input('周期与脉冲宽度之比M='); A=1;tao=T/M;a=tao/2;
y1=sym('Heaviside(t+a1)')*A; y=y1-sym('Heaviside(t-a1)')*A; y=subs(y,a1,a); y=simple(y);
%------------------------------ function x=time_fun_e(t)
% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 % t 是时间数组 % T 是周期 duty=tao/T=0.2 T=5;t=-2*T:0.01:2*T;tao=T/5;
x=rectpuls(t,1); %产生一个宽度tao=1的矩形脉冲 subplot(2,2,2) plot(t,x) hold on
x=rectpuls(t-5,1); %产生一个宽度tao=1的矩形脉,中心位置在t=5处 plot(t,x) hold on
x=rectpuls(t+5,1); %产生一个宽度tao=1的矩形脉,中心位置在t=-5处 plot(t,x)
title('周期为T=5，脉宽tao=1的矩形脉冲') axis([-10,10,0,1.2]) （d）运行结果及分析 ? 程序运行结果
调用CTFStpshsym.m函数文件，即可绘出周期矩形脉冲波形信号的双边幅度频谱。指令如下：在MATLAB命令窗口键入CTFStpshsym，回车，命令窗口将会出现：please Input所需展开的最高次谐波次数Nf=?；此处输入Nf=60，然后命令窗口将出现：please Input信号的周期T=?，此处输入T=5；然后命令窗口又出现：周期与脉冲宽度之比M=?，此处输入M=5，即可绘出频谱图，如图4-2所示。
29 图3-2 周期信号的双边频谱图