matlab基础学习记录之求傅里叶级数（stem的画图实现）

程序员文章站 2022-07-11 11:52:44

...

博主上次介绍了傅里叶级数，同时也介绍了三种基本函数的编写过程，本能的就想将傅里叶级数用matlab这个强大的工具来实现，但是在库里找了一圈也没有找到实现fourier seires的函数（或许是我没有认真找），但是抱着学习的心态，我还是自己试着去把傅里叶级数的函数实现写出来，这样，也能加深对matlab使用的理解。

万丈高楼平地起。

傅里叶级数，先把公式抛出来。

看起来是不是很简单？当初我也是这么想的，感觉三下五除二就能写出来，没啥难度。

但是事后发现，事情并没有想的那么简单。

下面，我就把过程中遇到的几个难点问题，以及解决方法，学习过程进行总结。

关于积分函数的选择

积分函数，通过查找资料，我发现主要有int、quad、quadl、integral，https://blog.csdn.net/u010999396/article/details/78861212在使用int做积分时，我发现，被积函数的原函数必须是存在的，int是先求原函数，再根据区间，求积分值。

例如在命令窗口对int进行试用。

当函数的原函数存在时

求函数的原函数（不定积分）

求该函数的定积分

下面用int求周期矩形脉冲（分段函数）的不定积分

报错！

而我们的被积函数，很明显，倘若f（t）是一个周期矩形脉冲，那么被积函数是不存在原函数的，因此，我们在求定积分时，就要采用数值积分，直接用类似求面积的原理来进行积分求解，将定积分问题，变成求和问题，避开了求原函数的过程。这里我们选用integral。

用int求周期矩形脉冲定积分（可以预想到肯定也报错）

下面我们使用integral函数，它不能求不定积分，却可以直接通过分解求和的方式求定积分。

因此在积分函数上，我们选择integral。

关于被积函数

积分函数的使用看起来很简单，但是我们发现，当被积函数是一个未知函数时，也需要我们认真考虑如何来构造积分函数。

使用integral函数求特定函数积分的例子在上面已经实践过，但是我们的被积函数是如下的形式,被积函数不在是一个特定的函数，而是由未知函数和已知函数相乘得到的新函数，那么，使用怎样的语句才能对这样一个新函数进行积分呢。

我想到的是，首先定义这个新函数，然后在命令行里运行积分函数，对新函数进行积分运算。

下面我们忘记傅里叶级数，使用一个简单的函数来验证我们的思路。

在这里未知函数使用为fun1 y=2t+p新函数结构为fun2

function y=fun1(t,p)
y=2.*t+p;
end

function y=fun2(n,f,p)
y=n+p.*integral(@(t)(f(t)+n),-p./2,p./2);
end

命令行

我们从积分函数的选择、被积函数的参数个数、被积函数的组成，都进行了简单函数的模拟，说明我们的思路是可行的，因此，下面将函数换成傅里叶级数

定义被积函数

function y=fs(n,f,p)    
%p为周期函数的周期
y=1./p.*integral(@(t)(f(t).*exp(-1i.*2.*pi.*n.*t./p)),-p./2,p./2);
end

但是在这里我又遇到了另外一个难点，在应用积分函数时，我们的被积函数y=fs(n,f,p)实际上有两个自变量n,t(p为函数周期会给出)，而积分时，除了被积的自变量，其他变量必须是已知的，因此，我们在上述命令行给n赋值为2，这样得到积分值18。

然而，为了画图，我们是想要的n未知的傅里叶级数，但如果想对f(t).*exp(-1i.*2.*pi.*n.*t./p)积分，就又必须给n赋值，这让我很矛盾。

由于被积函数里面有n，我根本无法通过integral函数，求出一个n仍为自变量的新函数，也就是我们无法利用这个新函数来stem画图。

这困扰了我很久，因此，我只能曲线救国了。

利用for循环求出多个n下的积分值，来进行stem画图。

脚本文件如下

clear all;
syms t n m;
p=10;
aaa@qq.com(t)zqjxmc(t,2,p);
for n=-20:1:20
    if abs(angle(fs(n,f,p)))>3.141 && angle(fs(n,f,p))<3.142
        m=-1;
    else
        m=1;
    end
    subplot(2,2,2);
    stem(n,abs(fs(n,f,p)));
    title('幅度谱');
    ylabel('|Fn|');
    xlabel('n(w=n*w0)');
    hold on;
    subplot(2,2,3);
    stem(n,angle(fs(n,f,p)));
    title('相位谱');
    ylabel('FAIn');
    xlabel('n(w=n*w0)');
    hold on;    
    subplot(2,2,4);
    stem(n,m.*abs(fs(n,f,p)));
    title('频谱');
    ylabel('Fn');
    xlabel('n(w=n*w0)');
    hold on;
end
q=-100:0.0001:100;
subplot(2,2,1);
plot(q,zqjxmc(q,2,p));
text(0,1.5,['脉冲宽度=2s']);
text(20,1.5,['周期=10s']);
title('周期矩形脉冲函数');
ylabel('f(t)');
xlabel('t(s)');
axis([-40,40,-1,2]);
grid on;

上面用到了一些画图的工具函数，之前都总结过。

占空比为1/5

matlab基础学习记录之求傅里叶级数（stem的画图实现）