【cg】【球谐光照】预备知识之傅立叶变换

前言

本篇会以多项式乘法作为引子，进一步了解快速傅立叶变换的原理和实现，从而深入理解傅立叶变换以及广义傅立叶变换的原理，为后面进一步推导和理解球谐函数打一下基础。

多项式乘法

已知

求

朴素的乘法

从

不难看出，如果我们把每一项展开两两相乘的话，拢共需要计算

次。

举个具体的例子。

设

则有

显然，总共进行了

次计算，最终得到的结果有

项，即最高次数为5次方。

这便是最朴素的多项式乘法了。接下来要做的事情就是让这件朴素的事情花哨起来，这就不得不从多项式的表示方式入手了。

多项式的表示

多项式可以用系数表示法和点值表示法来表示。

系数表示法 对于上面的多项式f(x)，可以使用一个系数向量来唯一确定，即

点值表示法 我们如果可以取f(x)上的n个自变量地方的值，那么这n个点值对也可以唯一确定这个多项式。

设

，则有

之所以这n个点值对可以唯一确定这个多项式，是因为我们将其分别代入f(x)，便可得到一个未知数为

的n元一次线性方程组，即

使用求解多元一次线性方程组的方法，如Gauss消元、LU分解、Crammer法则等[1]，便可求出该多项式所有的系数

。

通过这两种表示方式，我们可以寻找到求多项式乘法的另一种方式。

另一种乘法方式

现在我们已知的是

和

的系数表达式，即

我们需要求的是

的系数表达式，即

\[ F(x) \Leftrightarrow \vec{c} = (c_0, c_1, c_2, \ldots, c_{k-1}) \text{, (k = (n-1) + (m-1) + 1)} \]

其中

，即最高次项次数 + 1，即

的项数(包括系数为0的低次项)。

现在来考虑另一种思路，如果我们知道

上的

个点的值，组成

个点值对，那么我们便可以通过上面提到的解多元一次线性方程组，来得到

的系数表达式

。设这

个点值对为

很显然，对于任意定义域内的

，都有

那么

便可转化为

(注意这里的

并非向量积，而是向量的乘法，即各个分量相乘，得到相同项数的结果向量。)

那么我们只要分别在

和

上取

个点值对，组成它们的点值表示法，然后将它们相乘，便可得到被积函数

的点值表示法，然后再解多元一次线性方程组求出

的系数表达式。

还拿上面的

为例。

设

则要求的

为

\[ F(x) \Leftrightarrow \vec{c} = (c_0, c_1, c_2, \ldots, c_{k-1}) \text{, (k = (n-1) + (m-1) + 1) = 6} \]

那么我们接下来按上面的步骤来求一下

。

第一步 取

个点，即

，求得

和

的点值表达式，即

该步需要执行(

)次乘法，是因为通过秦九韶算法[2]，每一个

都只需要n次乘法得出，每一个

都只需要m次乘法得出。

根据系数表达式求其点值表达式的过程叫做求值。

第二步 根据上一步求出来的

和

的点值表达式，求出

的点值表达式，即

这一步只需要(

)次乘法。

第三步 将上同求出来的点值表示法，通过解多元一次线性方程组的方法转为系数表示法即可。

将

的元素分别代入

中，得

使用高斯消元法[3]转换成行阶梯式解方程组，可得

与上面朴素的乘法得到了相同的结果。

高斯消元的时间复杂度为

，由于在计算行阶梯式时每行的计算互不干扰，所以可使用多线程算法将复杂度降到

根据点值表达式求其系数表达式的过程叫做插值。

至此，我们通过求值(将系数转为点值)、点值相乘、插值(将点值转为系数)，三个步骤，消耗了比朴素的乘法

更高的复杂度

。

仿佛有些得不偿失-_-。

但是，这给优化时间复杂度提供了可操作的空间。因为如果我们选择一些特殊的点来进行插值和求值的话，便可将其复杂度降低到比朴素乘法更低的程度。

而单位复数根就是一种可以达成这种效果的特殊的点。

复数

定义

复数[4]是一个数域(C)，它是实数域(R)的延伸，额外包含带有虚数单位i的数。其中虚数单位i是-1的平方根，即

。

对于任意复数t都可记为

其中a、b都是实数，a称为t的实部，b称为t的虚部。所以，对于任意复数t,都可以由唯一有序对表示，即

在二维笛卡尔坐标系中，如果我们将x轴作用于实部a，将y轴作用于虚部b，便可在该坐标系下唯一确定一点与唯一一个复数一一对应，这个坐标系称为复平面。

我们也可以通过极角坐标系的方式表示复平面上的点。设

为与坐标系原点的距离，

为与x轴正方向的夹角且

，则有序对

也可唯一确定复平面上的一点，即可唯一确定一个与之对应的复数。如图所示。

复平面

复平面

即有

此时

即

同理也有

欧拉公式[5]将复数与三角函数关联起来，它的定义如下。

特别地，当

时，有欧拉恒等式

那么复数t可根据欧拉公式表示为指数形式，即

运算

设

。

加法

乘法

除此之外，由欧拉公式易得极坐标系下的乘法运算结果。

同理易得指数运算的结果如下。

单位复数根

单位复数是指

的复数，即落在复平面以坐标原点为圆心的单位圆上的复数。

单位复数根[6]的定义如下。

定义

n次单位复数根是指n次幂为1的复数，即

设

，则根据公式

有

故有

即对于给定的次数n,在

上共有n个单位复数根，它们是

记为

注意，如果k可以取到n，那么此时有

，它与

重复了，都是复平面上x轴正方向上的单位根，即

所以说n次单位复根总共只有n个，剩下的就是重复的了，这n个根在复平面上如图所示。

单位复数根

单位复数根

单位复数根动态

性质

单位复数根有一些后面计算需要用到的性质。

相消引理

由公式

可证。

折半引理

同样由公式

可得

在复平面上表现为旋转了180度，即方向相反的两个向量。

求和引理

可由等比数列求和公式证明。

在复平面上表现为所有的向量相加等于0。或者说所有向量终点组成的多边形的重心为坐标原点。

求和引理增强版

\[ \sum_{k=0}^{n-1} (\omega_n^k)^t = \left\{ \begin{array}{l} n & \text{, t \neq 0} \end{array} \right. \tag{3.6} \]

其中

。此公式也可由等比数列求和证明。

当

时，有

。

当

时，有

显然这是一个首项为1，公比为

的等比数列的前n项和。根据等比数列[9]的求和公式，可得

得证。

离散傅立叶变换DFT

离散傅立叶变换(DFT, Discrete Fourier Transform)[7]是指对于已知其系数表达式的n项多项式f(x)，分别代入这n个n次单位复数根

，得到n个点值对，形成该多项式的点值表达式，即

由

，若直接代入的话，时间复杂度为

。但是我们可以利用单位复数根的性质，将问题转化为规模更小的子问题，使用分治策略处理问题，可以使复杂度降至

，这也是使用单位复数根的原因。

我们首先来考虑一个项数为

的多项式f(x)，其中

为2的幂，即

。其实这已经考虑了所有的多项式，因为项数不是2的幂的多项式，总可以添加系数为0的项，一直加到项数为2的幂为止即可。

将其各项按照次数的奇偶性分组，即

提取奇次数项中的公因数

，可得两个关于

的多项式，分别记为

和

，即

\[ f(x) = \underbrace{(a_0 + a_2x^2 + \cdots + a_{n-2}x^{n-2})}_{\text{偶次数项A_1(x^2)}} + x \cdot \underbrace{(a_1 + a_3x^2 + \cdots + a_{n-1}x^{n-2})}_{\text{奇次数项A_2(x^2)}} \]

则有 \[ \left\{ \begin{array}{l} A_1(x) = (a_0 + a_2x + a_4x^2 + \cdots + a_{n-2}x^{\frac{n-2}{2}}),\text{共\frac{n}{2}项} \\ A_2(x) = (a_1 + a_3x + a_5x^2 + \cdots + a_{n-1}x^{\frac{n-2}{2}}),\text{共\frac{n}{2}项} \end{array} \right. \Rightarrow \begin{array}{l} f(x) = A_1(x^2) + x \cdot A_2(x^2) \tag{3.7} \end{array} \]

下面我们开始代入n个单位复数根

。

我们将k的取值区间分为左右两个部分，每个部分

个元素来分别讨论。

左侧部分 对于左边的

个数，即

，由公式

可得

由相消引理

可得

右侧部分 对于右边的

个数，我们将左边的

值向右偏移

，即仍然取

，通过代入

，取遍右边的值，即

。

此时由公式

有

由折半引理

可得

。

同时有

，再由相消引理

可得

。

所以上式可转化为

综上所述，我们的

值只需要取左侧的

个，即可通过

得到所有的

个点值对。

与此同时，我们只需要求出项数为

的多项式

的

个单位复数根

所表示的

个点值对即可，

同理。

从上面这句描述可以发现，这是一个与原问题完全相同，只不过规模从

变为

的一个子问题，所以我们完全可以使用分治的策略逐层划分为规则减半的子问题。

上述过程可以很容易转化成代码实现，一个有效的参考如下。

// complex
struct cp {
    double r, i;
    cp(double rr = 0, double ii = 0):r(rr), i(ii) {}
    cp operator + (cp x) { return cp(r+x.r, i+x.i); }
    cp operator - (cp x) { return cp(r-x.r, i-x.i); }
    cp operator * (cp x) { return cp(r*x.r-i*x.i, r*x.i+i*x.r); }
};

/*
* @input  a -- 系数组成的数组，长度为项数，a[i]表示第i项的系数。
* @input  n -- 项数。
* @output y -- 单位复数根处的值，y[i]表示第i个单位复数根处的值f(w_n^i)。
*/
void dft(in const cp a[], in int n, out cp y[]) {
    // 只有一个元素时为常数项，系数即为值
    if(n == 1) {
        y[0] = a[0];
        return;
    }

    // 第一步，系数分左右
    cp *ae = new cp[n>>1], *ao = new cp[n>>1];
    for(int i = 0; i < n/2; ++i) {
        ae[i] = a[i<<1];                        // A_1
        ao[i] = a[i<<1|1];                      // A_2
    }

    // 第二步，分治策略递归求规则更小的子问题
    cp *ye = new cp[n>>1], *yo = new cp[n>>1];
    dft(ae, n>>1, ye);                          // 求n/2个单位复根在A_1(x)的值
    dft(ao, n>>1, yo);                          // 求n/2个单位复根在A_2(x)的值

    // 第三步，合并结果，公式(3.10)
    cp wn = cp(cos(2*PI/n), sin(2*PI/n));       // 单位复数根的迭代增量
    cp w(1, 0);                                 // k = 0时的单位复根
    for(int i = 0; i < n/2; ++i) {
        y[i]       = ye[i] + w*yo[i];
        y[i + n/2] = ye[i] - w*yo[i];
        w = w*wn;                               // 相当于单位复根w_n^k中的k++
    }

    delete ae; delete ao; delete ye; delete yo;
}

显然它的时间复杂度为

，可通过主定理[8]求出。

离散傅立叶逆变换IDFT

离散傅立叶逆变换(IDFT, Inverse Discrete Fourier Transform)是指对于已知点值表达式的多项式f(x)，求其系数表达式

的过程。

IDFT同样可以使用单位复数根加速运算。

假设对于多项式

，

是我们需要求的系数表达式，已知其n个单位复数根所对应的点值对

，即

我们以这n个值

作为系数构造新的n项多项式g(x)，即

那么g(x)的表达式为

我们在g(x)上取n个点计算其值，这n个点满足如下条件，即上面单位复数根的共轭复数。

对于任意

，我们计算其值具体是多少。

代入

，得

设

，则可根据上面求和引理增强版

得

\[ G(j, k) = \left\{ \begin{array}{l} n & \text{, j \neq k} \end{array} \right. \]

那么上面的g(x)只有

即

项有值，即

那么对任意

，都有

所以我们只需要求出来

处的这n个g(x)的值，便是f(x)的系数表达式了，而求这n个g(x)的值的方法，就跟DFT一样使用分治策略即可，只不过把

就可以了，证明过程也与DFT一样，此处不作赘述。

总结一下流程。

1.已知f(x)在单位复数根

处的这n个值，将这n个值作为系数构造多项式g(x)。 2.求g(x)在

这n个地方的值，便可得到f(x)的系数表达式

。

这个过程与DFT只有两个差别，就是将

变成

，同时最终的结果除以

即可。所以稍微改一下DFT的代码即可实现。

/*
* @input  a -- 单位复数根处的值组成的数组，长度为项数，a[i]表示第w_n^i处的的值。
* @input  n -- 项数。
* @output y -- 系数*n，y[i]表示第i个系数*n,所以该结果要除以n才是真正的系数。
*/
void idft(in const cp a[], in int n, out cp y[]) {
    // 只有一个元素时为常数项，系数即为值
    if(n == 1) {
        y[0] = a[0];
        return;
    }

    // 第一步，系数分左右
    cp *ae = new cp[n>>1], *ao = new cp[n>>1];
    for(int i = 0; i < n/2; ++i) {
        ae[i] = a[i<<1];                        // A_1
        ao[i] = a[i<<1|1];                      // A_2
    }

    // 第二步，分治策略递归求规则更小的子问题
    cp *ye = new cp[n>>1], *yo = new cp[n>>1];
    idft(ae, n>>1, ye);                         // 求n/2个单位复根在A_1(x)的值
    idft(ao, n>>1, yo);                         // 求n/2个单位复根在A_2(x)的值

    // 第三步，合并结果，公式(3.10)
    cp wn = cp(cos(-2*PI/n), sin(-2*PI/n));     // 与dft唯一的修改 k -> -k
    cp w(1, 0);                                 // k = 0时的单位复根
    for(int i = 0; i < n/2; ++i) {
        y[i]       = ye[i] + w*yo[i];
        y[i + n/2] = ye[i] - w*yo[i];
        w = w*wn;                               // 相当于单位复根w_n^k中的k++
    }

    delete ae; delete ao; delete ye; delete yo;
}

可以使用额外的参数将DFT和IDFT整合起来。

/*
* @input  a -- 系数组成的数组，长度为项数，a[i]表示第i项的系数。
* @input  n -- 项数。
* @input  d -- 1表示DFT，-1表示IDFT。
* @output y -- DFT时表示单位复数根处的值，y[i]表示第i个单位复数根处的值f(w_n^i)。
* @output y -- IDFT时表示系数*n，y[i]表示第i个系数*n,所以该结果要除以n才是真正的系数。
*/
void fft(in const cp a[], in int n, out cp y[], int d = 1) {
    // 只有一个元素时为常数项，系数即为值
    if(n == 1) {
        y[0] = a[0];
        return;
    }

    // 第一步，系数分左右
    cp *ae = new cp[n>>1], *ao = new cp[n>>1];
    for(int i = 0; i < n/2; ++i) {
        ae[i] = a[i<<1];                        // A_1
        ao[i] = a[i<<1|1];                      // A_2
    }

    // 第二步，分治策略递归求规则更小的子问题
    cp *ye = new cp[n>>1], *yo = new cp[n>>1];
    fft(ae, n>>1, ye);                          // 求n/2个单位复根在A_1(x)的值
    fft(ao, n>>1, yo);                          // 求n/2个单位复根在A_2(x)的值

    // 第三步，合并结果，公式(3.10)
    cp wn = cp(cos(d*2*PI/n), sin(d*2*PI/n));   // 单位复数根的迭代增量
    cp w(1, 0);                                 // k = 0时的单位复根
    for(int i = 0; i < n/2; ++i) {
        y[i]       = ye[i] + w*yo[i];
        y[i + n/2] = ye[i] - w*yo[i];
        w = w*wn;                               // 相当于单位复根w_n^k中的k++
    }

    delete ae; delete ao; delete ye; delete yo;
}

迭代版实现

我们可以将上面的分治策略改为迭代实现。

注意到每次实施分治策略的时候，对于每一层递归，我们都将奇数项分在一起，偶数项分在一起。这相当于考察这一项在这一层递归中所对应的二进制位是否为1，如果是1，则说明它是奇数，分到奇数堆里面去，反正分到偶数堆里面去。

举个例子，当

时，对于一开始处于第三项的元素，它的二进制位是011，那么第一层递归的时候，考察它最右边的一位的值是1，所以它会被分到奇数堆里去，即下标会往右放，假设它最终被放到的下标的二进制位是xyz，那么x位必然为1,因为这一轮只有那4个偶数的x位才是0。

那么继续下一轮递归，此时考察它的中间的二进制位，这个二进制位也是1，那么它又被往右放了，且它最终被放到的下标的y位也是1。

继续下一轮，此时考察它的最左边的二进制位，这个二进制位是0，那么它也被往左放了，所以它最终下标的z位是0。

那么011到了最后一轮，便会被放到下标为110的位置，即它原来二进制位的反转结果。

如果我们确定了各个元素最终所在的下标，那么在每一层合并的时候都不会打乱上一层的下标顺序，就可以在原数组位置进行合并，如下代码所示。

// 二进制反转
int rev(int id, int n) {
    int ret = 0;
    for(int i = 0; (1<<i) < n; ++i) {
        ret <<= 1;
        if(id & (1<<i)) ret |= 1;
    }
    return ret;
}
// 另一种二进制反转
void rev(cp y[], int n) {
    for(int i = 1, j = (n>>1); i < n-1; ++i) {
        if(i < j) swap(y[i], y[j]); int k = (n>>1);
        while(j >= k) { j -= k; k >>= 1; } if(j < k) j += k;
    }
}
void FFT(cp a[], int n, cp y[], int d) {
    // 二进制反转 -- 先将元素放到正确的下标
    //for(int i = 0; i < n; ++i) y[i] = a[i]; rev(y, n); //另外一种二进制反转
    for(int i = 0; i < n; ++i) y[rev(i, n)] = a[i];

    // 从最底层开始迭代
    for(int i = 1; (1<<i) <= n; ++i) {                  // 层数
        int m = (1<<i);                                 // 当前层元素个数
        cp wn = cp(cos(d*2*PI/m), sin(d*2*PI/m));
        for(int j = 0; j < n; j += m) {                 // 每块的起始下标
            cp w = cp(1, 0);
            for(int k = 0; k < (m>>1); ++k) {           // 遍历每块的1/2个元素
                // 在原数组位置进行合并 -- 蝴蝶操作
                cp u = y[j+k];
                cp t = w*y[j+k + (m>>1)];
                y[j+k]          = u + t;
                y[j+k + (m>>1)] = u - t;
                w = w*wn;
            }
        }
    }

    // IDFT直接获得系数
    if(d == -1) for(int i = 0; i < n; ++i) y[i].r /= n, y[i].i /= n;
}

使用FFT实现多项式乘法

通过上面的fft，我们可以将插值和求值操作的复杂度都降到

，那么我们就可以使用fft来高效的实现多项式的乘法了。

// 返回大于x且离x最近的2的幂
int near(int x) {
    int i = 0;
    while(x > (1<<i)) ++i;
    return 1<<i;
}

inline int dcmp(long double a) { if(a < -eps) return -1; return (a > eps); }

// 多项式乘法
void poly_mul(double a[], int na, double b[], int nb, out double c[], out int& nc) {
    int i, j, n = na > nb ? na : nb;
    n = 1 << ((int)ceil(log(2*n)/log(2)-eps));  // n = near(na+nb-1);

    // f(x)和g(x)分别求值 -- O(nlogn)
    cp *x = new cp[n], *ya = new cp[n], *yb = new cp[n], *yc = new cp[n];
    for(i = 0; i < n; ++i) x[i] = cp(i < na? a[i] : 0, 0); FFT(x, n, ya, 1);
    for(i = 0; i < n; ++i) x[i] = cp(i < nb? b[i] : 0, 0); FFT(x, n, yb, 1);

    // 值相乘 -- O(n)
    for(i = 0; i < n; ++i) yc[i] = ya[i]*yb[i];

    // 插值 --得到最终结果 -- O(nlogn)
    FFT(yc, n, x, -1);

    for(i = 0; i < n; ++i) c[i] = x[i].r;
    for(nc = n; nc > 0 && dcmp(c[nc-1]) == 0; --nc);
    delete x; delete ya; delete yb; delete yc;

}

从代码可以看出，我们同样像文章开始那样通过求值(将系数转为点值)、点值相乘、插值(将点值转为系数)，三个步骤，但复杂度更低的

优化了多项式的乘法问题。这也是FFT最典型的应用之一了。

傅立叶级数

理解了离散傅立叶变换之后，我们可以继续更进一步理解更一般形式的傅立叶变换了。

接下来我们会先了解如何将周期函数通过傅立叶级数展开，然后再了解傅立叶级数是如何“进化”出傅立叶变换的，最后再了解又是怎么“究级进化”成广义傅立叶变换的。

在了解傅立叶级数之前，我们需要先了解一些信号统计和数学相关的前置知识。

周期信号和正弦波

周期信号广泛存在于我们的日常生活中，如简谐运动[10]、电磁波、引力波等。我们可以使用波动方程来描述周期信号对应的波动。其中最经典的要数正弦波[11]，其波动方程如下。

正弦波

正弦波

A是振幅，是振动的幅度。

是角频率(角速度)，若频率为

，则有

。

是相位，表示波在起始时间处的状态。 周期

我们可以线性组合多个不同参数的正弦波，以组成新的波动函数，从而形成新的波形。比如

这是一个自变量为时间的函数，显然，它会随着时间的偏移而波动，同时具有周期性，但是可能有不同于正弦波的波形，所以我们说它是一个新的波动函数。

值得注意的地方是，上面的波动函数，其实是一个时域上的函数，即它的自变量是时间，它的图像的x轴是以时间为单位的。

那么我们是不是可以以频率为自变量，或者说以角频率

为自变量，统计这个波在各个频率上的振幅。比如上面的函数，它在

的地方的振幅为1，在

的地方的振幅为2。那么当我们建立以

为x轴，以振幅为y轴的坐标系时，它的图像就变成了下图所示。

频域

频域

那么这个图像所表示的函数解析式便是一个离散的函数，即

\[ F(\omega) = \left\{ \begin{array}{l} 1 & \text{, \omega = 3} \\ \end{array} \right. \]

如果

的取值也是连续的，那么

就变成了一个连续的函数，这个函数的作用域是频率，所以这个函数我们又称之为频域上的函数。此时原来时域上的函数

也就变成了由无数个正弦波线性组合而成的了。

这就是喜闻乐见的时域与频域的概念。

频域与时域(图网侵删)

频域与时域(图网侵删)

除此之外，还有一个值得思考的问题，是不是所有的周期信号(周期函数)，都可以通过正弦波的线性组合来表示？或者更进一步，是不是所有的函数都可以通过正弦波的线性组合来表示？这个问题也许就是傅立叶当时思考的问题，我们后面再展开讨论。

级数

级数[12]是一个有穷或无穷序列的和。比如我们常见的等差数列和、等比数列和等。

其中笔者认为比较常用且需要掌握的便是泰勒级数。

泰勒级数

泰勒级数[13]是指在数学上，对于一个在实数或复数a邻域上，以实数作为变量或以复数作为变量的函数，并且是无穷可微的函数f(x)，它的泰勒级数是以下这种形式的幂级数

其中

表示

在点a处的n阶导数。如果

，则该级数又称麦克劳林级数。

如

的麦克劳林级数为

如果

在a处

次可导，那么对于这个区间的任意

，都有

该多项式称为

在a处的泰勒展开式，其中

是泰勒公式的余项，它是

的高阶无穷小。

可以看出，通过泰勒公式，我们可以将

以多项式的方式近似出来，这个多项式称为泰勒多项式。

下面我们开始研究傅立叶级数。

定义

对于满足一定条件的任意周期函数

，设其周期为

，角频率为

，则

都可以使用一系列的正弦函数展开，这里有一个比较炫酷的可视化效果网站[14]。

可以看出公式

中

为变量，

已知，振幅

和初相

是我们需要求的，但目前这种形式不太好求。

由三角恒等式两角和公式

，可得

同时设

(为了结果与

统一格式)，故有

这便是傅立叶级数的一种形式，最终求出的结果如下。

\[ \left\{ \begin{array}{l} a_n = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} f(t) \cdot \cos(n\omega t) dt & \text{, n \in \{1, 2, \ldots \}} \\ \end{array} \right. \tag{9.3} \]

公式

即为将

使用傅立叶级数展开的结果。

其中f(t)需要满足的条件如下。

f(t)在一个周期内连续，或在一个周期内只有有限个第一类间断点(左右极限都存在的间断点)。 f(t)至多只有有限个极值点。

上述条件又称狄利克雷充分条件。当满足上述条件时，我们说f(t)的傅立叶级数是收敛的，且有

若点t处f(t)是连续的，则级数收敛于f(t)。 若点t处f(t)是间断的，则级数收敛于

，其中

为无穷小。

下面来看上述结果是怎么推导出来的。

推导

观察公式

，我们能取到的所有的

函数值和

函数值，由以下集合组成。

这个集合里的元素在区间

正交，即两两乘积在这个区间的积分(专业术语叫做内积)为0，如下所示。

可通过三角恒等式证明，不作赘述。

我们称这个集合定义了一个坐标系三角函数系，集合里的元素都是该坐标系的原始轴(类比于迪卡尔坐标系的xyz)，我们称之为基或基底函数。那么我们求的

和

，就是f(t)在这个坐标系下的坐标，即f(t)在各个基底上的投影。

下面开始正式求解。

第一步 积分

设

，则由换元积分法,可得在

，即

上的积分为

所以有

所以有

第二步 乘

再积分

对公式

两边同乘以

，其中

，得

设

，则由换元积分法,可得在

，即

上的积分为

所以有

所以有

至此，我们求出了

的结果如上式所示，同时由第一步结果

也满足上式，故有

第三步 乘

再积分

第三步与第二步类似，只不过将

变成

即可。

对公式

两边同乘以

，其中

，得

设

，则由换元积分法,可得在

，即

上的积分为

所以有

所以有

至此，我们求出了

的结果如上式所示，即

至此，

得到了与一开始

一致的结果，这便是求解傅立叶级数的整个过程。

复数形式

根据欧拉公式

即

可得

代入

得

所以有

\[ \begin{array}{l} c_0 &=& \frac{a_0}{2} &=& \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) \mathrm{d}t \\ c_n &=& \frac{a_n-ib_n}{2} &=& \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) [\cos((-n)\omega t) +i\sin((-n)\omega t)] \mathrm{d}t = \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) e^{i(-n)\omega t} \mathrm{d}t & \text{, n = 1, 2, 3, \ldots} \\ \end{array} \]

显然这三个公式可以统一为以下形式。

\[ c_n = \frac{1}{T} \int\limits_{-\frac{T}{2}}^{\frac{T}{2}}f(t) e^{-in\omega t} \mathrm{d}t \text{, n = \ldots, -3, -2, -1, 0, 1, 2, 3, \ldots} \tag{9.8} \]

代入

，得

这便是傅立叶级数的复数形式了，十分优雅。

显然我们如果已知

，便可求出

。如果已知

，也可求出

，这也是时域和频域的相互转换，十分优雅。

不管是

中的

或

，还是

中的

，如果我们把它们都看作是自变量为角频率

，因变量为振幅的频域函数的话，那么每一个

都可以确认一个角频率

，显然这个频域函数的定义域显然是离散的，即有

其中

的取值是离散的。

所说如果f(x)的T逐渐增大的话，它的角频率

就会逐渐减小，那么相邻的

之间的距离就会减小，也就是说它的定义域的取值间隔越来越小，它的定义域在x轴上就会变得越来越密集。

考虑一种极端情况，当

时，

之间的间距

无穷小，此时相邻的

之间的距离为无穷小，即此时定义域将覆盖整个x轴，即

变成了自变量为频率的连续函数，且此时f(x)也不再是周期函数。而傅立叶变换考虑的就是这种情况。

傅立叶变换

定义

由上面的分析可以看出，傅立叶级数只对周期函数有效。而傅立叶变换[16]就是将傅立叶级数推广到非周期函数的方法。

观察复数形式的傅立叶级数

，我们上面也讨论了，当

时，

之间的间距

无穷小，所以我们从

得到了一个__连续__的频域函数，且

也并非周期函数了。

此时将有

证明

观察复数形式的傅立叶级数

，我们将

代入

，并取

，得

这里我们将

中的积分微元

换一个符号

，是为了与外层的

的自变量

作下区分。

显然，对于

积分项，外层的

是一个常数，可乘进积分项的每一项中，即

设

，即随着n的不同而在定义域坐标轴上选取的不同的位置。那么定义域上相邻两个自变量的距离为

代入

，得

这是一个

的极限，显然这等价于一个微元为

的积分，即

所以有

得证。

可以看到，离散的

变成连关于

的连续函数

。当我们传入一个频率作为自变量时，就会得到关于这个频率的信息(振幅、初相等)，所以又称

为时域函数

的频谱。可以直观的看出

的频率分布情况。

傅立叶变换将傅立叶级数扩展到非周期函数，这使我们具备了将任意非周期函数

从连续的时域函数得到连续的频域函数，或者从连续的频域函数得到连续的时域函数的能力。

但是对于周期函数，我们却无法使用傅立叶变换来展开，而是只能用傅立叶级数展开。

所以我们需找寻找一种可以统一傅立叶级数和傅立叶变换二者的公式，使傅立叶变换既可以应用于非周期函数，又可以应用于周期函数。这便是广义傅立叶变换。

广义傅立叶变换

在了解广义傅立叶变换之前，我们需要先来看一个特殊的函数，即狄拉克δ函数[17]，小名冲激函数。

冲激函数

定义 狄拉克δ函数的直观定义如下。

\[ \delta(x) = \left\{ \begin{array}{l} + \infty & \text{, x \neq 0} \\ \tag{11.1} \end{array} \right. \]

同时满足

筛选性质

函数有一个很重要的性质，即筛选性质即，即对于

定义域为

的函数

，有

根据定义

可知，当且仅当

时

项不为0，故

项可提出，得

之所以叫筛选性质，是因为它可以将

在

处的值筛选出来。

定义

对于

弦函数

，可以通过傅立叶变换得到其频域函数

如下。

同时频域函数

还可通过傅立叶逆变换变换回时域函数

，即

注意其中

为原函数的角频率，

为最终求得的频域函数的自变量。

与此同时，我们知道，对于任意周期函数，都可通过傅立叶级数展开成弦函数的线程组合，即公式

那么对于任意周期函数

，代入公式

，都可得到其频域函数

如下。

同理，对于任意周期函数的频域函数

，都可通过公式

对每一项进行傅立叶逆变换，也都可以重新得到其时域函数

，即公式

。

公式

即为广义傅立叶变换，它使得周期函数也可以通过傅立叶变换展开为

函数的线程组合。且只要证明

成立，即可推导出

。下面来证明。

证明

函数为一个非周期函数，我们对它进行傅立叶变换，即对时域函数

进行傅立叶变换，可得其频域函数

为

根据其筛选性质

可得

即

的频谱为一个常函数

。我们对该常函数重新进行傅立叶逆变换，便可得到原来的时域函数

，即

由于

和

的定义域都是

，所以对于上式，我们换个符号表示，即将

用

表示，将

用

表示，该式仍然成立，即

\[ \delta(\omega) = \frac{1}{2\pi}\int_{-\infty}^{\infty}e^{i\omega t} \mathrm{d}t \text{, \omega \in (-\infty, \infty)} \tag{11.7} \]

由

取

显然公式

成立。

由

函数的筛选性质

，易证公式

也成立，即

得证。

后记

至此，我们了解了离散傅立叶变换、傅立叶级数、傅立叶变换以及广义傅立叶变换的一些并不十分严谨的定义和证明。正如前文所述，本文的侧重点是对概念和思想的熟悉，以及推导过程中的重要步骤和结果，顺带着科普一些基础但很有趣的数学知识，所以有选择地忽略了很多细枝末节(如边界条件、级数的敛散性及相关证明、一些性质和应用等)，读者可自行脑(goo)补(gle)。

个人认为最重要的是了解时域和频域概念及其互相转换的思想。我们工作中也许会碰到一些在时域中处理不了或者不好处理的问题，转换到频域上处理可能能找到一条更好的道路。

我们在傅立叶级数的推导过程中，有提到过内积的概念，以及三角函数系的概念。其实这里的三角函数系在泛函分析中是一个希尔伯特空间，即完备的内积空间，它相当于将三维的欧几里德空间推广到n维乃至无限维但又不失完备性。而傅立叶变换相于在n维的希尔伯特空间的上向各个坐标轴(基底)的投影。

而广义傅立叶变换与泛函数分析就更密切了，因为

函数本身就不是一个严格意义上的函数。它可以定义为一个分布，即广义函数，即给定一个

，

可以得到某个值，即

是测试函数

的一个线性泛函。

函数可以使用测度的概念严谨定义，它其实是n维希尔伯特空间的一个测度，一个函数相对于

函数的积分便可认为是相对于这个测度的勒贝格积分。

以上两段没有什么实际参考意义，纯粹是为了显得专业而胡说八道，限于笔者自身水平，本文的绝大部分公式的推导和证明都是初等的微积分而已，应该很容易理解。不过如果读者掌握过《测度论》《泛函分析》等更专业的知识，也许就可以从另外一个更加高级的抽象层次理解和描述问题，无奈本人非数学专业出身，也没读过研-_-||，目前来说很难再从更高一层的视角作抽象和总结，希望以后可以找个时间填一下这些坑，此刻只能深表遗憾了。

参考

[1]多元线性方程组的解

[2]秦九韶算法

[3]高斯消元法

[4]复数

[5]欧拉公式

[6]单位复数根

[7]离散傅立叶变换

[8]主定理

[9]等比数列

[10]简谐运动

[11]正弦波

[12]级数

[13]泰勒级数

[14]Fourier-Visualizations

[15]An Interactive Introduction to Fourier Transforms

[16]傅立叶变换

[17]狄拉克δ函数