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); !次谐波叠加 f21=trifourierseries(a0, an, bn, 21, t); E次谐波叠加 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()中即可可以跑得比较快，不过这样的话在修改为不同函数进行积分时，这里也要进行修改，各有所长，也有所缺，看个人所好以及实际情况进行取舍即可。