浅尝辄止傅里叶

目录

• 前言
• 从音乐到声波
• 信号
• 傅里叶级数
• 傅里叶变换
  • 离散傅里叶变换
  • 卷积定理
  • 多项式卷积
  • 快速傅里叶变换
  • 快速数论变换
    • 时域抽取(DIT)
    • 频域抽取(DIF)
  • 转置原理
    • 多项式多点求值
• 从傅里叶级数到画图
• 音频/图片存储中的傅里叶变换
• 后记

前言

不会复变函数，不会傅里叶分析，不会信号与系统，没进省队没有国赛牌，大概是写得不太好，如题目所言：浅尝辄止，数学不行，不写证明，语文太弱，通顺就算赢。
抽象的概念比较多，因为有些时候没有必要特意展开说。
如有硬伤请务必指出。

图中许多的周转圆绘制出了 \(\aleph\)，下面将探讨这是怎么做到的。

从音乐到声波

随着自然科学的进步，人们以机械波这个模型描述声音，即常用的“声波”的概念。于是有几个描述声波性质的名词：音调（频率），响度（振幅）以及音色。
不同的乐器音色不同，但是只要音域有交集就可以演奏出相同音调的单音。而多个单音可以组成和弦，直观感受类似单音，于是可以发现，声波在这里叠加了。如何把叠加的声音分离？
下面引入一个抽象模型来更好的解决这些问题。

信号

信号（一维信号）可以被理解为一个关于时间的函数，描述的时候也有许多与函数类似的描述，例如周期性信号就是这个关系具有周期性。对于上面提到的声波，我们可以将其看作一个时间-气压的函数。

连续时间信号：在时间上连续。
离散时间信号：相当于取数个点值。

正弦型信号是一种常见的信号，形式为 \(f(t) = A \sin(\omega t + \varphi)\)，其中 \(A\) 为振幅，\(w\) 为角频率，\(\varphi\) 为初相位。正弦型信号是周期信号，周期为 \(T = \dfrac{2\pi}{\omega} = \dfrac{1}{f}\)，\(f\) 为频率。
根据欧拉公式\(e^{ix} = \cos x + i \sin x\)，可得：

\[\large{ sin(\omega t) = \dfrac{1}{2i}(e^{i\omega t} - e^{-i\omega t})\\ cos(\omega t) = \dfrac{1}{2}\ (e^{i\omega t} + e^{-i\omega t}) } \]

傅里叶级数

在为热传导进行数学建模的过程中，傅里叶提出，一切周期信号（函数）都可以用许多的不同初相，振幅，频率的正弦波的叠加来表示，如下式：

\[\large f(t) = A_0 + \sum_{n = 1}^{\infty} A_n \sin(n\omega t+\varphi_n) \]

这里可以发现每个正弦信号的角频率都是 \(\omega\) 的整数倍。
于是我们得到了这个级数，只需将其逐层求出，我们就能得到用正弦信号的叠加表示之，傅里叶通过麦克劳林级数转化上式，结论如下：

\[\large \begin{aligned} & f(t) = \dfrac{a_0}{2} + \sum_{n = 1}^{\infty} [a_n\cos(n\omega t) + b_n\sin(n\omega t)]\\ & a_n=\dfrac{2}{T}\int_{t_0}^{t_0 + T} f(t) \cos(n \omega t)\\ & b_n=\dfrac{2}{T}\int_{t_0}^{t_0 + T} f(t) \sin(n \omega t) \end{aligned} \]

能够以如此表示一个函数，要点在于三角函数的正交性：

\[\large{ \int_{-\pi}^{\pi} \cos(nx)\ \mathrm{d}x = 0\ (n \in \mathbb{N^{*}})\\ \int_{-\pi}^{\pi} \sin(nx)\ \mathrm{d}x = 0\ (n \in \mathbb{N^{*}})\\ \int_{-\pi}^{\pi} \cos(nx)\cos(kx)\ \mathrm{d}x = 0\ (n,k \in \mathbb{N^{*}},n\ne k)\\ \int_{-\pi}^{\pi} \cos(nx)\sin(kx)\ \mathrm{d}x = 0\ (n,k \in \mathbb{N^{*}},n\ne k)\\ \int_{-\pi}^{\pi} \sin(nx)\sin(kx)\ \mathrm{d}x = 0\ (n,k \in \mathbb{N^{*}},n\ne k) } \]

一组正交的向量是线性无关的，将其作为一组基底就可以表出其张成的空间中所有的向量。

从上文我们可以得到一个特性，举个例子，我们有一个信号，要求其中有没有频率为 \(\omega\) 的部分，我们可以将这个信号与频率为 \(\omega\) 的正弦信号相乘，然后求一个周期内的积分，如果积分的值较大，即频率为 \(\omega\) 的信号占原来信号的一部分，因为根据正交性，如果频率不相同的，相乘后积分，一个周期内得到的就是 \(0\)。

傅里叶变换

有一句话，叫做：“傅里叶变换是时域到频域的变换”，围绕这句话我们做一点展开。
根据傅里叶级数，信号可被正弦波的组合表出，若取一个极小的 \(\omega\)，其整数倍可以视作是连续的，于是可以求出任何频率的正弦信号在一个信号中分离出多大“程度”。

原函数（信号）是定义在时间上的，通过一维的时间而确定其意义，这是时域，而现在根据上述方法定义一个新函数，其自变量是频率 \(\omega\)，得到的是这个频率“占用”的多少，这个新函数就是频域上的，通过原函数得到这一特定新函数的过程，称之为傅里叶变换。

下图中蓝色图像为信号时域的图像，后方的三个是分离出的正弦信号的时域图像，以 \(f_0,2f_0\) 为轴的为该信号频域的图像。

举个例子，乐谱上的音符描述的是音调（频率），而演奏出的是声波，其中乐谱是频域信息，演奏的声波波形是时域信息，无论演奏出的声波在波峰还是波谷，乐谱都是一样的，频域和时域不互相干扰。

为了完成这一过程，根据上面傅里叶级数的一节提到的，将其与特定频率的函数相乘，然后求积分：

\[\large{ F(\omega) = \mathcal{F}[f(t)] = \int_{-\infty}^{\infty} f(t) e^{-i\omega t}\ \mathrm{d}t\\ f(t) = \mathcal{F}^{-1}[F(\omega)] = \dfrac{1}{2\pi}\int F(\omega) e^{i\omega t}\ \mathrm{d}t }\]

对于一个非周期函数，只需要将其当作周期为 \(+\infty\) 看待，所以上式是对任意函数的傅里叶变换。

3B1B 采用了一种特殊的理解方式，通过“缠绕”来理解。
下面用文字复述一下：
观察式子中的 \(e^{-i\omega t}\)，相当于在复平面绕一圈后，恰好能按照给定的频率把函数的两个周期的图像“缠绕”在单位圆上（即三角函数的复数形式，复平面上转动，加负号相当于顺时针）
在缠绕之后，如果这个频率不合适，那么原函数的尖峰会比较均匀地分布在圆周上，现在我们取其中一些点，这些点的质心越接近原点，说明这个函数与这一次取的频率越不契合（参见正交性那部分）。

于是这个远离原点的程度就是能从原函数里分离出这个频率的正弦波的程度。

\[\huge\begin{aligned} &\color{Purple}{F}(\color{Green}{f}) = \color{Plum}{\sum_{n=0}^{N-1}} \color{RoyalBlue}{f(n)} \color{OrangeRed}{e}^{\color{OrangeRed}{-i} \color{Orange}{2\pi} \color{Green}{f} \color{Plum}{\frac{n}{N}}}\\ &欲求\color{Green}{某一频率}\color{Purple}{能从原函数分离出多少}，\\ &先\color{Green}{按照那个频率}\color{OrangeRed}{旋转}\color{RoyalBlue}{原函数}\color{Orange}{一周}，\\ &再求出\color{Plum}{所有点贡献的和}。 \end{aligned}\]

然后我们将采样数量上升到无穷多，即使用积分，就得到了傅里叶变换。
由此可知，傅里叶变换具有线性性。直观地理解，先将信号叠加后分解，和先分解后叠加，效果相同。

离散傅里叶变换

对于一些信号，虽然是连续的，但是我们不能找出一个合适的数学表达式去描述，这时候我们可以选择取样，取出足够多的数据点也可以代表那个信号。
对于一个序列 \(\{x_n\}\)，其离散傅里叶变换为

\[\large X_{f} = \sum_{n = 0}^{N-1} x_n e^{-i2\pi f \frac{n}{N}} \]

同样，还有逆变换

\[\large x_{n} = \dfrac{1}{N}\sum_{n = 0}^{N-1} X_n e^{i2\pi f \frac{n}{N}} \]

离散傅里叶变换在时域与频域都离散。

卷积定理

离散傅里叶变换可以用于计算卷积，因为其有一条性质如下：

\[\large{ f_1(x) \leftrightarrow F_1(\omega),f_2(x) \leftrightarrow F_2(\omega)\\ \mathcal{F}[f_1(x) * f_2(x)] = F_1(\omega) \cdot F_2(\omega) }\]

也就是可以将两个函数分别傅里叶变换后求其频域乘积，然后逆变换得到的函数就是原先两函数的卷积。
证明：

\[\large\begin{aligned} \mathcal{F}[f_1(x) * f_2(x)] &= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f_1(\tau)f_2(x - \tau)\ \mathrm{d} \tau \ e^{-i\omega x}\ \mathrm{d} x\\ &=\int_{-\infty}^{\infty} f_1(\tau) \int_{-\infty}^{\infty} e^{-i\omega (x - \tau)} f_2(x - \tau)\ \mathrm{d}xe^{-i\omega \tau} \ \mathrm{d} \tau\\ &=F_2(\omega)\int_{-\infty}^{\infty}f_1(\tau) e^{-i \omega \tau}\ \mathrm{d}\tau\\ &= F_1(\omega)*F_2(\omega) \end{aligned}\]

多项式卷积

用 \(n + 1\) 个不同点的点值可以表出唯一的 \(n\) 次多项式，于是可知，多项式的傅里叶变换可以进行 \(n + 1\) 个采样然后采用离散傅里叶变换。
对于 \(n\) 次多项式 \(A(x)\) 与 \(m\) 次多项式 \(B(x)\)，其卷积为一个 \(n + m - 1\) 次多项式，于是只须求出 \(n + m\) 个点值即可对应唯一的 \(A(x) * B(x)\)。

然后我们回到离散傅里叶变换：

\[\large X_{f} = \sum_{n = 0}^{N-1} x_n e^{-i2\pi f \frac{n}{N}} \]

这里的 \(e^{-i2\pi f \frac{n}{N}}\) 是旋转一周然后按频率进行均分。
引入一个概念：复单位根。

考虑如下方程 \(x^n = 1 (n \in \mathbb{N}^*)\)，只考虑实数域的解时，\(x = 1\) 显然成立，在 \(n\) 为偶数的时候 \(x = -1\) 也是一个根。现在把范围放到复数域，可以得到复数根形式如下：\(e^{\frac{2\pi ki}{n}}(k=\mathbb{N}^*)\)，记为 \(\omega_{n}^{k}\)，放到复平面上，就是把单位圆均分为 \(n\) 份，用直线连接就得到了一个复平面上的正 \(n\) 边形。

复单位根有以下性质：

• \(\omega^{k}_{n} = \omega^{k + n}_{n}\)
• \(\omega^{k}_{n} = \omega^{tk}_{tn} (t \in \mathbb{N}^*)\)
• \(\omega^{k + n}_{2n} = -\omega^{k}_{2n}\)

那么如果我们对一个多项式 $A(x) = a_{n - 1}x^{n - 1} + a_{n - 2}x^{n - 2} + \cdots + a_0 $ 的系数序列做傅里叶变换，得到的就是 \(A(e^{-\frac{2\pi ik}{n}})\) 这个点值，因此傅里叶变换可以快速求出多项式代入复单位根得到的一系列特殊点值。

然后考虑逆变换的过程，发现区别只在于一个 \(\frac{1}{N}\) 以及 \(e\) 上面的指数把负号摘掉了：

\[\large x_{n} = \dfrac{1}{N}\sum_{n = 0}^{N-1} X_n e^{i2\pi f \frac{n}{N}} \]

这时候取单位根的倒数，最后乘上 \(\frac{1}{N}\) 即可还原到多项式的系数形式。
上面提到了傅里叶变换是一个线性变换，于是离散傅里叶变换可以被描述为一个线性算子，以矩阵表示：

\[\large \begin{aligned} \begin{bmatrix} 1 &1 &1 &\cdots &1\\ 1 &\omega_{n}^{1} &\omega_{n}^{2} &\cdots &\omega_{n}^{n - 1}\\ 1 &\omega_{n}^{2} &\omega_{n}^{4} &\cdots &\omega_{n}^{2(n - 1)}\\ \vdots &\vdots &\vdots &\ddots &\vdots\\ 1 &\omega_{n}^{n - 1} &\omega_{n}^{2(n - 1)} &\cdots &\omega_{n}^{(n - 1)(n - 1)} \end{bmatrix} \begin{bmatrix} a_0\\a_1\\a_2\\ \vdots \\a_{n - 1} \end{bmatrix} = \begin{bmatrix} v_0\\v_1\\v_2\\ \vdots \\v_{n - 1} \end{bmatrix} \end{aligned} \]

逆变换同理。

下面给出一种不需要卷积定理的“傅里叶变换能做卷积”的解释。
用点值表示多项式，称之为点值表示法，对于两个点值表示的多项式，我们只需要 \(\mathcal{O}(n)\) 地把对应的点值相乘就能得到结束后的点值，于是如果我们能求出点值和通过点值复原多项式就可以求其卷积。
直接计算 \(n\) 个点值和通过拉格朗日插值复原多项式是 \(\mathcal{O}(n^2)\) 的，复杂度瓶颈就在这里，但是现在我们可以通过离散傅里叶变换求一个多项式的特殊点点值表示，但是按照 DFT 公式来计算仍然是 \(\mathcal{O}(n^2)\)的，于是现在需要一种低于 \(\mathcal{O}(n^2)\) 的方法来进行离散傅里叶变换，就可以求多项式卷积，这一方法就是经常被提到的 快速傅里叶变换。

（本文是在已经从傅里叶级数得到连续以及离散傅里叶变换的公式的基础上推出的 DFT 的理解方式，只从单位复根角度也可以理解 DFT 的过程，可以参考：https://www.cnblogs.com/RabbitHu/p/FFT.html ）

快速傅里叶变换

FFT 的基本思想是分治法，也就是对于一个一次的多项式，可以 \(\mathcal{O}(1)\) 地求出一点值，配合一定的拆分，合并的手段，将复杂度降低到 \(\mathcal{O} (n\log n)\)

快速数论变换

将运算放在一个模大质数的域中，此时一个原根就有类似复单位根的性质，用原根代替复单位根，乘法逆元代替取倒数，就完成了 FFT 向 FNTT 的转变，优点是不需要进行浮点运算，理解 DFT 就能理解 NTT，不赘述。

时域抽取(DIT)

如何把对整个多项式的 FFT 分为更小规模的问题？首先要能分成两个相似的形式，为了每次都能恰好对半分开，我们可以在高次不断补 \(0\) 使系数个数为 \(2\) 的幂。
观察一个多项式 \(A(x) = a_{n - 1}x^{n - 1} + a_{n - 2}x^{n - 2} + \cdots + a_0\)，按照时域上的奇偶，即多项式次数分成两组，一侧指数为 \(1,3,5,7,\cdots\)，一侧指数为 \(0,2,4,6,\cdots\)，可以发现对应位置相差 \(1\)，只需要把奇数部分提取一个 \(x\)，两部分就相同了，分解为两个小规模的问题：

\[\large\begin{aligned} X_{k} &= \sum_{n = 0}^{N - 1} x_n e^{-i2\pi k \frac{n}{N}} (k = 0,1,2,\cdots,N - 1)\\ &= \sum_{n 为偶数} x_n \omega_{N}^{nk} + \sum_{n 为奇数} x_n \omega_{N}^{nk}\\ &= \sum_{m = 0}^{\frac{N}{2} - 1} x_{2m} \omega_{N}^{2mk} + \omega_{N}^{k}\sum_{m = 0}^{\frac{N}{2} - 1} x_{2m + 1} \omega_{N}^{2mk}\\ &= \sum_{m = 0}^{\frac{N}{2} - 1} x_{2m} \omega_{N/2}^{mk} + \omega_{N}^{k}\sum_{m = 0}^{\frac{N}{2} - 1} x_{2m + 1} \omega_{N/2}^{mk}\\ \end{aligned}\]

而 \(\large\sum_{m = 0}^{\frac{N}{2} - 1} x_{2m}\omega_{N/2}^{mk}\) 和 \(\large\sum_{m = 0}^{\frac{N}{2} - 1} x_{2m + 1} \omega_{N/2}^{mk}\) 就是分离出两部分的 DFT 的结果，简记为：

\[\large X_{k} = X_1[k] + \omega_{N}^{k}X_2[k] \]

上面推出的是对于频域前 \(\dfrac{N}{2}\) 求出的，现在考虑后半：

\[\large\begin{aligned} X_{k + \frac{N}{2}} &= X_1[k + \frac{N}{2}] + \omega_{N}^{k + \frac{N}{2}}X_2[k + \frac{N}{2}]\\ &= X_1[k] - \omega_{N}^{k}X_2[k] \end{aligned}\]

于是我们得到了 DIT-FFT。

考虑对其进行优化，在分为子问题的过程中，需要按照奇偶分组，那么在进行 FFT 之前把每一个系数放到变换完成后其位置，就可以直接向上合并，跳过了分隔的过程。
发现在 FFT 后每个数到了自身原先下标二进制按位翻转（总位数为 \(\lceil\log_2(N)\rceil - 1\)，不足在前面补 \(0\)）后的位置，二进制翻转可以 \(\mathcal{O}(n)\) 递推得到，形如：
rev[i] = ((rev[i >> 1] >> 1) | ((i & 1) << (n - 1)));
这个操作名为蝴蝶变换/位逆序变换。

频域抽取(DIF)

现在我们按照频域的奇偶进行分组：
偶数部分：

\[\large\begin{aligned} X_{2k} &= \sum_{n = 0}^{N - 1} x_{n} \omega_{N}^{2k\cdot n}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} x_n \omega^{2k\cdot n}_{N} + \sum_{n = \frac{N}{2}}^{N-1} x_{n} \omega^{2k\cdot n}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} x_n \omega^{2k\cdot n}_{N} + \sum_{n = 0}^{\frac{N}{2}-1} x_{n + \frac{N}{2}} \omega^{2k\cdot (n + \frac{N}{2})}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} x_n \omega^{2k\cdot n}_{N} + \sum_{n = 0}^{\frac{N}{2}-1} x_{n + \frac{N}{2}} \omega^{2k\cdot n}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} (x_n + x_{n + \frac{N}{2}}) \omega^{kn}_{\frac{N}{2}} \end{aligned}\]

奇数部分：

\[\large\begin{aligned} X_{2k + 1} &= \sum_{n = 0}^{N - 1} x_{n} \omega_{N}^{(2k + 1)\cdot n}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} x_n \omega^{(2k + 1)\cdot n}_{N} + \sum_{n = \frac{N}{2}}^{N-1} x_{n} \omega^{(2k + 1)\cdot n}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} x_n \omega^{(2k + 1)\cdot n}_{N} + \sum_{n = 0}^{\frac{N}{2}-1} x_{n + \frac{N}{2}} \omega^{(2k + 1)\cdot (n + \frac{N}{2})}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} x_n \omega^{(2k + 1)\cdot n}_{N} - \sum_{n = 0}^{\frac{N}{2}-1} x_{n + \frac{N}{2}} \omega^{(2k + 1)\cdot n}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} (x_n - x_{n + \frac{N}{2}}) \omega^{(2k + 1)n}_{N}\\ &= \sum_{n = 0}^{\frac{N}{2}-1} \omega_{N}^{n}(x_n - x_{n + \frac{N}{2}}) \omega^{kn}_{\frac{N}{2}} \end{aligned}\]

令：

\[\large\begin{cases} x_1[n] = x_n + x_{n + \frac{N}{2}}\\ x_2[n] = \omega_{N}^{n}[x_n - x_{n + \frac{N}{2}}] \end{cases}\]

则有：

\[\large\begin{aligned} &X_{2k} = \sum_{n = 0}^{\frac{N}{2}-1}x_1[n]\omega_{\frac{N}{2}}^{kn}\\ &X_{2k + 1} = \sum_{n = 0}^{\frac{N}{2}-1}x_2[n]\omega_{\frac{N}{2}}^{kn} \end{aligned}\]

同样也可以分为两个子问题。

由上图所示，DIT-FFT 将蝴蝶变换后的系数转化为点值，而 DIF-FFT 将系数转化为蝴蝶变换后的点值，二者互为转置，因此在进行 DFT 的时候使用 DIF-FFT，在进行 IDFT 的时候使用 DIT-FFT 就可以省去位逆序变换，在预处理 \(\omega\) 的时候也按照蝴蝶变化后的顺序处理，能在此基础上减少移动 \(\omega\) 数组指针的次数。
然后我们就得到了一种比较快的多项式卷积实现方式，即使用 DIF-DIT 省略蝴蝶变换，用 FNTT 代替 FFT，用迭代代替递归。

\(\Large 代码在这里\)

转置原理

(大概需要一点线性代数基础？)
首先定义线性算法为如下过程：
有一 \(n\) 阶方阵 \(A\)，由一个 \(n\) 行的列向量 \(a\)，求出 \(b = Aa\)。
为了方便表示可以将矩阵乘法展开：

\[\large b_i = \sum_{i,j} A_{i,j}a_j \]

对于其转置，即形如 \(b' = A^T a'\) 的过程。
展开后就是

\[\large b_i =\sum_{i,j}A_{j,i}a_j \]

转置原理：任何线性算法都有转置算法，复杂度不变。
于是优化原算法可以转化为优化转置后的算法。

然后考虑一般化转化为转置算法的过程：

1. 单位矩阵 \(I\) 可以通过初等变换得到任意可逆矩阵。
2. 一次初等变换可以改写为乘一个初等矩阵（\(I\) 经过一次初等变换得到的矩阵）

不同的初等矩阵效果不同：

1. 将 \(I\) 的第 \(i\) 行与第 \(j\) 行交换，乘一个列向量可以将其第 \(i\) 与第 \(j\) 位交换。
2. 将 \(I\) 的 \((i,i)\) 位置的值从 \(1\) 改为 \(a\)，乘一个列向量可以使其第 \(i\) 位乘 \(a\)。
3. 将 \(I\) 的 \((i,j)\ (i \ne j)\) 位置的值从 \(0\) 改为 \(a\)，乘一个列向量可以使其第 \(j\) 位乘 \(a\) 的值加到第 \(i\) 位上。

于是我们尝试将 \(A\) 表示为 \(A = E_1E_2\dots E_k\)
然后转置：\(A^T = E_k^T E_{k - 1}^T\dots E_1^T\)
也就是进行所有初等变换的顺序反序，前两个操作不影响，第三种变为使其第 \(i\) 位乘 \(a\) 的值加到第 \(j\) 位上。

前文提到 DFT 可以写成矩阵形式，该矩阵对称，其转置就是自身，以及 DIT-FFT 和 DIF-FFT 互为转置。

现在观察两个多项式 \(a,b\) 的卷积：

\[\large a * b = \sum_{i = 0}^{n + m - 2} \sum_{j = 0}^{i} a_i b_{j - i} x^i \]

次数 \(n\) 的多项式与次数 \(m\) 的多项式运算得到次数 \(n + m - 1\) 的多项式。
这里将 \(b\) 当作常量的话，卷积过程可以视为一个线性算法，对于其转置，就是给出 \(n - 1\) 次多项式和 \(m - 1\) 次多项式，得到一个 \(n - m\) 次多项式。
形如：

\[\large (a*b)^T = \sum_{i = 0}^{n + m - 2} \sum_{j = \max(0,i - n + m)}^{\min(i,m - 1)} f_{i} g_{j} x^{i - j} \]

然后这个东西是 \(a\) 与系数反转后的 \(b\) 卷积后的第 \(m - 1\) 到 \(n - 1\) 位，因为溢出的部分并不需要，所以转置后常数减小。

多项式多点求值

范德蒙德矩阵为一个 \(n \times n\) 的方阵，形式如下：

\[\large{ V_a = \begin{bmatrix} 1 &a_0 &a_0^2 &\cdots &a_0^{n - 1}\\ 1 &a_1 &a_1^2 &\cdots &a_1^{n - 1}\\ 1 &a_2 &a_2^2 &\cdots &a_2^{n - 1}\\ \vdots &\vdots &\vdots &\ddots &\vdots\\ 1 &a_{n - 1} &a_{n - 1}^2 &\cdots &a_{n - 1}^{n - 1}\\ \end{bmatrix}}\]

如果将 \(a_0,a_1,\cdots,a_{n- 1}\) 视为变量序列，那么对于一个 \(n - 1\) 次多项式 \(b\)，\(V_a b\) 就是在 \(a_0,a_1,\cdots,a_{n- 1}\) 这些点得到的多项式点值。

为了优化多点求值的过程，取 \(V_a\) 的转置：

\[\large V_a^T = \begin{bmatrix} 1 &1 &1 &\cdots &1\\ a_0 &a_1 &a_2 &\cdots &a_{n - 1}\\ a_0^2 &a_1^2 &a_2^2 &\cdots &a_{n - 1}^2\\ \vdots &\vdots &\vdots &\ddots &\vdots\\ a_0^{n - 1} &a_1^{n - 1} &a_2^{n - 1} &\cdots &a_{n - 1}^{n - 1}\\ \end{bmatrix}\]

那么转置后求的就是 \(v_i = \sum_{j = 0}^{n - 1} a_{j}^i b_j\)。
如果将 \(v_i\) 视作一个多项式的系数，得到：

\[\large\begin{aligned} V(x) &= \sum_{i = 0}^{n - 1} (\sum_{j = 0}^{n - 1} a_j^i b_j) x^i \\ &= \sum_{i = 0}^{n - 1} \sum_{j = 0}^{n - 1} (a_j^i b_j) x^i \\ &= \sum_{j = 0}^{n - 1} b_j \sum_{i = 0}^{n - 1} a_j^i x^i \end{aligned}\]

然后后面是等比数列求和，即：

\[\large V(x) = \sum_{j = 0}^{n - 1} \frac{b_j}{1 - a_j x} \]

于是我们需要计算这个多项式前 \(n\) 项的系数，将和式通分，得到：

\[\large \dfrac{\sum_{i = 0}^{n - 1} b_i \prod_{j \ne i}(1 - a_j x) }{\prod_{k = 0}^{n - 1}(1 - a_k x)} \]

令 \(F_{l,r}\) 表示项数 \([l,r)\) 区间内的分子，\(G_{l,r}\) 表示项数 \([l,r)\) 区间内的分母，有：

\[\large{ G_{l,r} = G_{l,mid}G_{mid,r}\\ F_{l,r} = F_{l,mid}G_{mid,r} + F_{mid,r}G_{l,mid} }\]

分治计算 \(F_{0,n},G_{0,n}\)，最后一步多项式求逆即可。

目前我们得到的是多项式多点求职的转置的解法，现在考虑将上述过程再转置。

1. 从下到上求转为从上到下求。
2. 加法卷积转置，变为减法卷积。

此时根节点的 \(F'\) 就是输入的多项式和 \(\dfrac{1}{G_{0,n}}\) 减法卷积。
因为转置后到了一个分治位置 \([l,r)\) 需要立即使用 \(G_{l,r}\) 这时候我们需要先处理出 \(G\)，然后用 \(F'_{l,r}\) 更新两个子区间：

\[\large{ F'_{l,mid} = (F'_{l,r}G_{mid,r})^T\\ F'_{mid,r} = (F'_{l,r}G_{l,mid})^T }\]

递归到底，就得到了点值。

\(\Large 代码在这里\)

突然意识到如果不用 vector 存储就会MLE，好麻烦啊
还有一个由 FFT 衍生的 Bluestein 算法，但是笔者太菜了并不会，可以尝试阅读 2020 年集训队论文。

从傅里叶级数到画图

回看周转圆绘图，如果我们在时域上分析，那么绘图过程可以视作一个周期函数 \(f(t)\)，对于输入一个时间 \(t\) 有唯一的复平面上的点，也就是 \(f(t) = g(t) + i\cdot h(t)\)，在 \(t\) 时刻点的坐标就是 \((g(t),h(t))\)。

然后我们看一个函数：\(e^{ix}\)，

可以发现以某个频率旋转的点的运动，在实轴上投影为 时间-实部 的余弦函数，在虚轴上投影为 时间-虚部 的正弦函数（可以从欧拉公式里看出来），现在回看傅里叶级数，每一项展开，有其频率，初相，将每一项视为圆周运动，其运动的叠加就构成了被绘制的图形，实际计算过程中，计算机使用的是离散傅里叶变换，在图像上取点，然后得到每个圆周运动的频率。

音频/图片存储中的傅里叶变换

MP3 格式在存储的过程中，会使用傅里叶变换将分解出对整体影响更小的小幅度高频部分，不存储这部分信息，以压缩文件大小。

JPEG格式图像在存储的时候，首先把图像分解成许多 \(8\times8\) 大小的像素块，RGB 用三个通道分别表示。可以发现，如果一块像素块颜色非常接近，其频域的分量将大部分为 \(0\) ，只有出现了的颜色的部分有值，于是可以通过一轮傅里叶变换（这里其实是离散余弦变换，傅里叶变换衍生出的）就可以减少大量相同信息的占用。

同样的，对声音/图像的频域分析还可以用于滤波，这里也使用了傅里叶变换。（以及很多衍生出的变换，比如拉普拉斯变换，小波变换）

后记

就是说如果有人对我说：“我研究一条一头加热了的金属棒子研究出一个能分解几乎所有周期函数的通法，还能解微分方程”，我肯定是不信的。
啥都能缝的神经网络已经缝上 DFT 了，用神经网络做 DFT,神秘。

感谢我自己，克服了拖延症坚持查资料把这篇破玩意写完了。
感谢 fxj，虽然他在写这篇博客方面没有帮我，但是我就是想提一句，膜拜大佬。

傅里叶级数画图的动图使用 3B1B 的数学动画引擎 Manim 制作。
时域-频域图片来自网络，yandex以图搜图找到的比较明晰的版本。
\(e^{i2\pi ft}\) 图像图片来自 Understanding Digital Signal Processing's Frequency Domain
DIT-DIF 图像来自维基百科。
DIT-DIF 优化的 FFT 代码参考了 yhx-12243 与 クトリ 的多项式模板。
转置原理参考 2020 集训队论文 《转置原理的简单介绍》 陈宇。