画傅立叶级数

假设函数 $s(x)$ 定义在 $[x_0, x_0+P]$ 的区间，前面提到 $s(x)$ 可写成其傅立叶级数的和：

$$s(x) = \frac{a_0}{2} + \sum_{n=1}^{\infty} (a_n\cos(\frac{2\pi n x}{P}) + b_n\sin(\frac{2\pi n x}{P})),$$

其中

$$a_n = \frac{2}{P}\int_{x_0}^{x_0+P} s(x)\cos(\frac{2\pi n x}{P})\ dx$$

$$b_n = \frac{2}{P}\int_{x_0}^{x_0+P} s(x)\sin(\frac{2\pi n x}{P})\ dx$$

$a_n$ 及 $b_n$ 为 $n$ 阶谐波的系数。

接下来，我们要使用 Octave 来观察一个锯齿状的三角波，我们把它的傅立叶级数的系数算出来，然后把 $N$ 项谐波的成分加起来，和原来波形放在一起观察。 这边使用的波形 $s(x)$ 函数如下：

可以算出其傅立叶级数为：

$$s(x) = \frac{2}{\pi}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\sin(nx)$$

程式步骤如下：

1. MATLAB 有一个 sawtooth 函数，可以直接画出锯齿波形。sawtooth 定义在 Signal Processing 工具箱，如果是使用 MATLAB，必须确定此工具箱已正确安装；如果使用 Octave，则可以使用 signal 套件。关于 Octave 套件的安装方式，在 Windows 环境下，可参考 官网的说明，如果是在 Ubuntu 的话，直接使用以下指令进行安装：

   sudo apt-get install octave-signal

2. 在 Octave 里面，先使用 pkg load signal 指令将 signal 套件载入。接著便可以使用 sawtooth 函数（可用 help sawtooth 查看此函数的使用说明）。试著使用以下指令画出 sawtooth 函数图形：

   t = -8:0.1:8;
   s = sawtooth(t-pi); % sawtooth 预设周期为 2*pi，在 0 的地方转折
   plot(t,s);
   grid on; % 画格线

   观察一下波形看是否正确。

3. 画部分级数和。

   N = 5;   % 项数
   sgn = 1; % 用作正负变号使用，第一项为正
   t = -8:0.1:8;
   fs = zeros(size(t)); % 设 fs 跟 t 维度一样，全为 0
   for n = 1:N          % 从 1 跑到 N
       fs = fs + 2*sgn*sin(n*t)/(n*pi); % 加一项
       sgn = sgn * -1;  % 变号
   end
   plot(t,fs); % 画图

   观察一下合成的波形。

4. 最后将两个波形画在一起，程式码如下：

   N = 5;
   t = -8:0.1:8;
   sgn = 1;
   s = sawtooth(t-pi); % sawtooth 预设周期为 2 pi，在 0 的地方转折
   fs = zeros(size(s));
   for n = 1:N
       fs = fs + 2*sgn*sin(n*t)/(n*pi);
       sgn = sgn * -1;
   end
   plot(t,s,'b',t,fs,'r');
   grid on; % 画格线

   执行后的结果如下：

   

在上面的程式中，我们已经知道 $s(x)$ 的傅立叶级数的 $a_n$ 和 $b_n$ 的系数，如果是其他函数的话，可能要花很多时间计算。实际上，Octave 有一个 Symbolic Math Toolbox，可以帮忙做符号运算。

如果是使用 Octave，可以先用以下指令安装 symbolic 套件：

sudo apt-get install octave-symbolic

接著用 pkg load symbolic 载入套件来使用。

例如 $f(x)=x^2$，假设定义在 $[-\pi,\pi]$ 的区间，可以使用以下方式计算 $a_n$ 的系数。

pkg load symbolic
syms x n;    % x, n 为符号
f = x^2;     % f(x) = x^2
an = int(f*cos(n*x), x, -pi, pi) / pi;    % 计算定积分的值

練習 3

Octave 里面有一个 square 函数，可以用来画出方波函数 (用 help square 查询看看)。请计算方波的傅立叶级数 (也可上网查询或使用 symbolic 套件计算)，接著用程式画出方波及其部分傅立叶级数和的图形，此处方波的周期和振幅可自己定义。