Matlab求解周期函数的傅⾥叶级数以及作频谱图与相位图（⼀）前⾔
Matlab并没有⾃带的求解傅⾥叶级数的函数，本⽂将介绍如何使⽤Matlab进周期函数的傅⾥叶级数分析，内容包括：
1、求解傅⾥叶级数的系数
2、求N次谐波的叠加函数，画图⽐较与原函数的差值
3、做出傅⾥叶级数的幅度谱与相位谱
（⼆）傅⾥叶级数系数的求解
设f(x)为周期为T的周期函数，则我们有傅⾥叶级数展开式：
根据系数的求解的定义，我们使⽤int()函数进⾏积分即可求解，如果f(x)在⼀个周期内为分段函数的话可能还需分段积分，我是⾃⼰写了⼀个统⼀处理，当然也可以⾃⼰分段写，本质是⼀样的。这⾥以⼀个周期三⾓函数为例进⾏求解，三⾓波函数图像如下：
则在⼀个周期内的函数表达式为：
则求解系数的代码如下：
syms t n;
T=2;w=2*pi/T;
f1=2*t+1;f2=-2*t+1;
a0=1/T*myint('t',f1,-0.5,0,f2,0,0.5);
an=myint('t',f1*cos(n*pi*t),-0.5,0,f2*cos(n*pi*t),0,0.5);
bn=myint('t',f1*sin(n*pi*t),-0.5,0,f2*sin(n*pi*t),0,0.5);
运⾏结果为：
（三）作N次谐波的合成图并与原函数进⾏⽐较
作合成图实际上就对多个函数进⾏⼀定项数的叠加，为了适应不同项数的叠加处理，这⾥编写函数进⾏实现：
%trifourierseries.m
function [ f ] = trifourierseries( a0, an, bn, m, t )
%TRIFOURIERSERIES 求傅⾥叶级数m次谐波的合成
%  a0、an、bn为傅⾥叶级数的系数
%  t为变量(取样间隔也就是⾃变量)
f=a0;
syms n;
for n=1:m
f=f+eval(an)*cos(n*pi.*t);
f=f+eval(bn)*sin(n*pi.*t);
end
接着就是进⾏画图处理：
%求傅⾥叶变换
t=-6:0.01:6;
d=-6:2:6;
%tripuls为三⾓波函数，进⾏偏移叠加处理可以得到⼀个类似三⾓周期函数的图
f=pulstran(t,d,'tripuls');
%3次谐波叠加
f3=trifourierseries(a0, an, bn, 3, t);
%9次谐波叠加
f9=trifourierseries(a0, an, bn, 9, t);
%21次谐波叠加
f21=trifourierseries(a0, an, bn, 21, t);
%45次谐波叠加
f45=trifourierseries(a0, an, bn, 45, t);
%级数展开图
subplot(2,3,1);plot(t,f,'r',t,f3,'b');grid on
axis([-6.1,6.1,-0.1,1.1]);title('展开3项');
xlabel('t');ylabel('f(t)');
subplot(2,3,4);plot(t,f,'r',t,f9,'b');grid on
axis([-6.1,6.1,-0.1,1.1]);title('展开9项');
xlabel('t');ylabel('f(t)');
subplot(2,3,2);plot(t,f,'r',t,f21,'b');grid on
axis([-6.1,6.1,-0.1,1.1]);title('展开21项');
xlabel('t');ylabel('f(t)');
subplot(2,3,5);plot(t,f,'r',t,f45,'b');grid on
axis([-6,6,-0.1,1.1]);title('展开45项');
xlabel('t');ylabel('f(t)');
（四）幅度谱与相位图作图
先给定需要绘制的范围，再对具体的幅度An以及相位ψn进⾏求解，最后画出An~n以及ψn~n即可：
%幅度谱--相位谱
n=0:1:10;
anVal=eval(an);
bnVal=eval(bn);
%注意A0需要⾃⼰赋值
An=sqrt(anVal.^2+bnVal.^2);An(1)=a0;
%phi0同理
phi=atan(-bnVal./anVal);phi(1)=0;
subplot(2,3,3);stem(n,An,'b');
grid on; axis([-0.1,10.1,-0.1,1.1]);
title('幅度谱');xlabel('n');ylabel('An');
subplot(2,3,6);plot(n,phi,'b');
grid on; axis([-0.1,10.1,-0.1,1.1]);matlab求傅里叶变换
title('相位谱');xlabel('n');ylabel('ψn');
（五）最终绘图结果
（六）说明
该例⼦为周期三⾓波函数的求解，如需改成其它函数，只需要将最开始的分段f1、f2以及相应的积分过程进⾏修改即可。
该程序直接运⾏的话会⽐较慢，因为在进⾏傅⾥叶级数的计算是将a0、an、bn进⾏传参然后求解，既涉及到符号运算⼜涉及到了数值运算，故运算特别慢，在算出a0、an、bn以后，直接将三者的值代⼊trifourierseries()中即可可以跑得⽐较快，不过这样的话在修改为不同函数进⾏积分时，这⾥也要进⾏修改，各有所长，也有所缺，看个⼈所好以及实际情况进⾏取舍即可。