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]);
title('相位谱');xlabel('n');ylabel('ψn');

（五）最终绘图结果

（六）说明

该例子为周期三角波函数的求解，如需改成其它函数，只需要将最开始的分段f1、f2以及相应的积分过程进行修改即可。

该程序直接运行的话会比较慢，因为在进行傅里叶级数的计算是将a0、an、bn进行传参然后求解，既涉及到符号运算又涉及到了数值运算，故运算特别慢，在算出a0、an、bn以后，直接将三者的值代入trifourierseries()中即可可以跑得比较快，不过这样的话在修改为不同函数进行积分时，这里也要进行修改，各有所长，也有所缺，看个人所好以及实际情况进行取舍即可。