第 30 章 多项式与快速傅里叶变换

本章对应原书第 30 章。量化研究里到处是"卷积"：移动平均是收益序列和一个权重核的卷积，自相关函数是序列和自己的相关，独立损失之和的分布是各自分布的卷积。直接算长度 \(n\) 的卷积要 \(\Theta(n^2)\) 次乘法；快速傅里叶变换（FFT）把它降到 \(\Theta(n\lg n)\)。原书用"多项式乘法"这个干净的模型把 FFT 讲清楚：多项式的系数相乘就是卷积，而 FFT 是在系数表示与点值表示之间快速切换的工具。本章先按原书把算法和证明讲透，再在量化实战里把它用到滚动统计、周期检测、损失分布和期权定价上。

学习目标

读完本章，你应当能够：

1. 说清多项式的系数表示与点值表示各自擅长什么运算，并解释"求值—逐点相乘—插值"三步为什么能把乘法从 \(\Theta(n^2)\) 降到 \(\Theta(n\lg n)\)。
2. 证明单位复根的消去引理、折半引理、求和引理，并说明每一条在 FFT 和逆 FFT 中起什么作用。
3. 写出递归版和迭代版 FFT（蝴蝶操作、位逆序置换），推导 \(T(n)=2T(n/2)+\Theta(n)\)，并能用 FFT 计算逆 DFT。
4. 正确使用卷积定理：知道什么时候必须补零，什么时候得到的是循环卷积，以及 numpy 的符号约定与原书有何不同。
5. 在量化场景中用 FFT 计算任意核的滚动加权统计、自相关函数、周期图和组合损失分布，并识别"全样本频域滤波"带来的前视偏差。
6. 了解 FFT 期权定价（Carr–Madan）的结构：一次 FFT 得到整条执行价网格上的价格。

读前导读

这一章在解决什么问题

这一章只解决一个计算问题：卷积怎么算得快。"卷积"听着陌生，其实你天天在用：

• 移动平均、EWMA、任何"对过去 \(L\) 期收益加权求和"的指标，都是收益序列和一组权重的卷积。
• 信用组合里，每笔贷款的损失是一个小的离散分布；独立贷款加总后的损失分布，就是这些小分布一个接一个卷积出来的。你在 CFA 里算"两笔贷款各自违约或不违约，共有四种情况"，就是一次最小的卷积。
• 掷两颗骰子求点数和的分布，也是卷积。

直接算卷积，长度 \(n\) 的两个序列要做大约 \(n^2\) 次乘法。\(n=1000\) 时是一百万次，\(n=10^6\) 时是一万亿次，这就慢到不能用了。快速傅里叶变换（FFT）把它降到大约 \(n\log_2n\) 次：\(n=10^6\) 时约两千万次，快了几万倍。

原书讲 FFT 的角度是"多项式乘法"。原因是：把序列的各项当作多项式的系数，两个多项式相乘时，乘积的系数恰好就是两个序列的卷积。于是"快速算卷积"变成"快速乘多项式"。快速乘的办法分三步：先把两个多项式在一批特殊的点上求值，再把对应的值逐个相乘（这一步很便宜），最后从乘积的值反推出系数。FFT 的全部巧妙之处，在于选了一批特殊的点（单位复根），让"求值"和"反推"都能用分治法快速完成。

你不需要会推导 FFT 才能用它——numpy.fft 和 scipy.signal.fftconvolve 一行就能调用。本章对你真正有用的是三件事：知道它为什么快（从而知道什么问题可以用它加速）；知道必须补零，否则序列尾部会"绕回"开头，在回测中造成前视偏差；知道"全样本频域滤波"不能当作实时可得的信号。

需要先想起来的数学

1. 复数与欧拉公式。 复数写成 \(a+bi\)，其中 \(i\) 满足 \(i^2=-1\)。可以把它看作平面上的点 \((a,b)\)。乘法规则是 \((a+bi)(c+di)=(ac-bd)+(ad+bc)i\)。关键直观：模长为 1 的复数可以写成 \(e^{iu}=\cos u+i\sin u\)（欧拉公式），它是单位圆上角度为 \(u\) 的点；乘以 \(e^{iu}\) 就是把一个点绕原点转 \(u\) 角。例如 \(i=e^{i\pi/2}\)，乘以 \(i\) 就是转 90°，所以 \(i\cdot i\) 转了 180°，等于 \(-1\)。欧拉公式可以由 \(e^x\)、\(\cos x\)、\(\sin x\) 的幂级数展开得到，见 第 00 册第 04 章 级数与收敛。

2. 等比数列求和。 \(1+q+q^2+\cdots+q^{n-1}=\frac{q^n-1}{q-1}\)（\(q\ne1\)）。你在年金现值公式里用过它。本章的"求和引理"就是把 \(q\) 换成一个复数后套这个公式。见 第 00 册第 04 章。

3. 矩阵乘向量与矩阵的逆。 在 \(n\) 个点上求一个多项式的值，可以写成"一个 \(n\times n\) 矩阵乘系数向量"；反推系数就是乘这个矩阵的逆。只要知道"矩阵可逆 ⟺ 行列式不为零 ⟺ 方程组有唯一解"即可。见 第 00 册第 06 章 线性代数速成。

4. 分治递归式 \(T(n)=2T(n/2)+\Theta(n)\)。 意思是"把规模 \(n\) 的问题拆成两个规模 \(n/2\) 的子问题，再花与 \(n\) 成正比的时间合并"。拆 \(\log_2n\) 层、每层总工作量约 \(n\)，合计约 \(n\log_2n\)。本册 第 04b 章 递归式求解与主定理 有完整讨论；归并排序是同一个式子。（\(\Theta(\cdot)\) 表示"同阶"，\(\lg\) 表示以 2 为底的对数。）

怎么读这一章

核心必读：30.1（卷积是什么）、30.2.4（"求值—相乘—插值"的整体方案）、30.5.3（补零与循环卷积）、30.8.1 和 30.8.4（滚动统计和前视偏差）。这几节读懂了，就能在实际工作中正确使用 FFT。

建议顺序：先读 30.1 和 30.2，理解整体方案；然后跳到 30.5.3 和 30.8 的量化实战，看它实际怎么用；最后有兴趣再回头读 30.3–30.4 的单位复根和 FFT 分治，那里是"为什么快"的数学原因。

第一次可以只看结论的部分：30.3 三条引理的证明（记住"单位复根平方后个数减半""单位复根之和为 0"两句话即可）；30.6 的高效实现（位逆序、迭代版、并行电路是工程细节）；30.7 的思考题选讲；30.8.5 Carr–Madan 期权定价（需要特征函数的背景，可在读完第 08 册后回来看）。

30.1 从多项式乘法说起

30.1.1 多项式与卷积

域 \(F\)（本章通常取复数域 \(\mathbb C\)）上的多项式（polynomial）是形式和

两个次数界为 \(n\) 的多项式相加，只要把系数逐项相加：\(c_j=a_j+b_j\)，\(\Theta(n)\) 时间。原书的例子：

相乘则麻烦得多。乘积 \(C(x)=A(x)B(x)\) 的次数界为 \(2n-1\)：

式 (30.2) 定义的系数向量 \(c\) 称为 \(a\) 与 \(b\) 的卷积（convolution），记作 \(c=a\otimes b\)。一般地，\(\deg C=\deg A+\deg B\)；次数界分别为 \(n_a\)、\(n_b\) 的多项式之积次数界为 \(n_a+n_b-1\)，原书为了方便常直接说 \(n_a+n_b\)。

推导拆解：用最小的例子看清式 (30.2)。取 \(A(x)=1+2x\)，\(B(x)=3+x\)，系数向量 \(a=(1,2)\)、\(b=(3,1)\)。展开 \((1+2x)(3+x)=3+x+6x+2x^2=3+7x+2x^2\)。

按 (30.2) 逐项核对：\(c_0=a_0b_0=3\)；\(c_1=a_0b_1+a_1b_0=1\cdot1+2\cdot3=7\)；\(c_2=a_1b_1=2\)。规律是：\(c_j\) 收集所有"下标加起来等于 \(j\)"的乘积 \(a_kb_{j-k}\)，因为 \(x^k\cdot x^{j-k}=x^j\)。

把系数换成概率就是分布的卷积：若 \(a=(0.9,0.1)\) 是"贷款 1 损失 0 或 1 单位"的概率，\(b=(0.8,0.2)\) 是"贷款 2 损失 0 或 1 单位"的概率，那么 \(c=(0.72,\,0.26,\,0.02)\) 就是两笔独立贷款总损失为 0、1、2 的概率。\(c_1=0.9\times0.2+0.1\times0.8\)，正是"恰好一笔违约"的两种情形之和。

30.1.2 卷积在量化里的样子

把式 (30.2) 换一个角度看。设 \(r_0,r_1,\dots\) 是收益序列，\(w_0,\dots,w_{L-1}\) 是一组权重，那么

本章所有记号中 \(i\) 专指虚数单位 \(\sqrt{-1}\)，不作下标用。

30.2 多项式的两种表示

30.2.1 系数表示

系数表示（coefficient representation）就是系数向量 \(a=(a_0,a_1,\dots,a_{n-1})\)（视为列向量）。它擅长两件事：

• 求值：用 Horner 法则（Horner's rule）在 \(\Theta(n)\) 时间内算出 \(A(x_0)\)：
  \[A(x_0)=a_0+x_0\Big(a_1+x_0\big(a_2+\cdots+x_0(a_{n-2}+x_0a_{n-1})\cdots\big)\Big).\]
• 加法：\(\Theta(n)\)。

def horner(a, x0):          # a[0..n-1] 为升幂系数
    y = 0
    for j in range(len(a) - 1, -1, -1):
        y = a[j] + x0 * y
    return y                # 时间 Θ(n)，空间 O(1)

但乘法要 \(\Theta(n^2)\)。

30.2.2 点值表示

次数界为 \(n\) 的多项式也可以用 \(n\) 个点值对来表示：

从系数得到点值，就是在 \(n\) 个点上求值（evaluation），逐点用 Horner 法则需 \(\Theta(n^2)\)；后面会看到，点选得巧可以降到 \(\Theta(n\lg n)\)。反过来，从点值恢复系数叫插值（interpolation）。插值是否总能做、结果是否唯一？

定理 30.1（插值多项式的唯一性） 对任意 \(n\) 个点值对 \(\{(x_k,y_k)\}\)，只要 \(x_k\) 两两不同，就存在唯一的次数界为 \(n\) 的多项式 \(A\) 满足 \(A(x_k)=y_k\)，\(k=0,\dots,n-1\)。

证明：条件 \(y_k=A(x_k)\) 写成矩阵方程

用 LU 分解解这个方程组要 \(O(n^3)\)。更快的是拉格朗日公式（Lagrange's formula）：

白话解释：定理 30.1 是中学几何"两点确定一条直线"的推广：两个点确定一条直线（次数界 2），三个点确定一条抛物线（次数界 3），\(n\) 个横坐标不同的点确定唯一一个次数界为 \(n\) 的多项式。所以"写出全部系数"和"写出在 \(n\) 个点上的值"是同一个多项式的两种等价描述，可以互相转换。

这个想法你在收益率曲线上见过：既可以用参数（如 Nelson–Siegel 的几个系数）描述一条曲线，也可以用若干关键期限上的收益率描述它。前者便于在任意期限上求值，后者便于逐点比较和加减。多项式的两种表示也是这样各有所长，下一小节说明为什么点值表示做乘法特别便宜。

"病态"（ill-conditioned）的意思是：输入的微小误差会被放大成输出的巨大误差。例如用很多等距点拟合高次多项式，个别点的数据稍有扰动，系数就可能剧烈变化。

30.2.3 两种表示下的运算

点值表示擅长的运算：

• 加法：若 \(A\)、\(B\) 用同一组点表示，\(C=A+B\) 的点值就是 \((x_k,y_k+y'_k)\)，\(\Theta(n)\)。
• 乘法：\(C=AB\) 的点值就是 \((x_k,y_ky'_k)\)，也只要 \(\Theta(n)\)。

乘法有一个细节：\(C\) 的次数界是 \(2n\)，而唯一确定一个次数界 \(2n\) 的多项式需要 \(2n\) 个点（习题 30.1-4：少于 \(n\) 个点不能唯一确定次数界 \(n\) 的多项式）。所以 \(A\)、\(B\) 必须先用 \(2n\) 个点表示，称为扩展点值表示（extended point-value representation）。

点值表示的弱项是在新的点上求值——没有比先转回系数表示再用 Horner 更简单的办法。

两种表示各有长短，于是自然的想法是：乘法时先转到点值表示，逐点相乘，再转回来。关键在转换能否快。

30.2.4 快速乘法的整体方案

原书图 30.1 给出的方案如下。设 \(n\) 是 2 的幂（不是就补高位零系数）。

1. 倍增次数界：给 \(A\)、\(B\) 各补 \(n\) 个高位零系数，变成次数界 \(2n\) 的多项式。\(\Theta(n)\)。
2. 求值：各做一次长度 \(2n\) 的 FFT，得到它们在 \(2n\) 个 \(2n\) 次单位复根处的值。\(\Theta(n\lg n)\)。
3. 逐点相乘：得到 \(C\) 在这 \(2n\) 个点上的值。\(\Theta(n)\)。
4. 插值：对这 \(2n\) 个点值做一次逆 DFT（同样用 FFT），得到 \(C\) 的系数。\(\Theta(n\lg n)\)。

定理 30.2 两个次数界为 \(n\) 的多项式可以在 \(\Theta(n\lg n)\) 时间内相乘，输入和输出都是系数表示。

整个方案的成败在第 2、4 步：在任意 \(2n\) 个点上求值要 \(\Theta(n^2)\)，插值更贵。秘诀是选单位复根作求值点——在那里求值就是离散傅里叶变换（DFT），插值就是逆 DFT，二者都能用 FFT 在 \(\Theta(n\lg n)\) 内完成。下面两节把这件事讲透。

30.3 单位复根

30.3.1 定义

\(n\) 次单位复根（complex \(n\)th root of unity）是满足 \(\omega^n=1\) 的复数 \(\omega\)。恰好有 \(n\) 个：

称为主 \(n\) 次单位根（principal \(n\)th root of unity），其余单位根都是它的幂。原书脚注说明：很多信号处理文献取 \(\omega_n=e^{-2\pi i/n}\)，数学实质相同——numpy、scipy 正是用负号约定，30.6 节会再提醒。

\(n\) 个单位根 \(\omega_n^0,\dots,\omega_n^{n-1}\) 在乘法下构成一个群，结构与加法群 \((\mathbb Z_n,+)\) 相同：\(\omega_n^j\omega_n^k=\omega_n^{(j+k)\bmod n}\)，\(\omega_n^{-1}=\omega_n^{n-1}\)。可以把它想成一块表盘：乘 \(\omega_n\) 就是指针转动 \(1/n\) 圈。

白话解释：以 \(n=4\) 为例把单位复根具体写出来。\(\omega_4=e^{2\pi i/4}=e^{i\pi/2}=\cos90°+i\sin90°=i\)。它的各次幂是 \(\omega_4^0=1\)，\(\omega_4^1=i\)，\(\omega_4^2=i^2=-1\)，\(\omega_4^3=-i\)，\(\omega_4^4=1\)（转满一圈回到起点）。这四个数在复平面上恰好是单位圆上的东、北、西、南四个点。

为什么偏偏选这些点来求值？它们有两条普通实数点没有的性质。第一，平方后个数减半：\(1,i,-1,-i\) 平方后是 \(1,-1,1,-1\)，只剩 \(\{1,-1\}\) 两个值，这让问题能一分为二（折半引理）。第二，加起来抵消：\(1+i+(-1)+(-i)=0\)，这让逆变换的公式非常简单（求和引理）。若改用 \(1,2,3,4\) 这样的实数点，平方后是 \(1,4,9,16\)，个数不减，分治就做不下去。

如果你没学过复数，可以先把 \(\omega_n^k\) 当成"一个满足 \((\omega_n^k)^n=1\) 的特殊数"，只记住上面两条性质，就能跟上后面的推导。

30.3.2 三条关键性质

引理 30.3（消去引理，cancellation lemma） 对任意整数 \(n\ge0,k\ge0,d>0\)，

推论 30.4 对偶数 \(n>0\)，\(\omega_n^{n/2}=\omega_2=-1\)。（转半圈就是 \(-1\)。）

引理 30.5（折半引理，halving lemma） 若 \(n>0\) 为偶数，则 \(n\) 个 \(n\) 次单位复根的平方恰好是 \(n/2\) 个 \(n/2\) 次单位复根，每个出现两次。

证明：由消去引理，\((\omega_n^k)^2=\omega_n^{2k}=\omega_{n/2}^k\)。再看 \(k\) 与 \(k+n/2\) 两个根：

折半引理是分治的核心：它保证"对一半规模的子问题求值"时，求值点的个数也正好减半。

引理 30.6（求和引理，summation lemma） 对任意整数 \(n\ge1\) 和不被 \(n\) 整除的非零整数 \(k\)，

直观上，把 \(n\) 个等距分布在圆上的向量相加，它们互相抵消。求和引理将用于证明逆 DFT 公式。

30.4 DFT 与 FFT

30.4.1 DFT 的定义

我们要在 \(n\) 个 \(n\) 次单位复根 \(\omega_n^0,\dots,\omega_n^{n-1}\) 上对次数界为 \(n\) 的多项式

例（原书习题 30.2-2） 求 \(a=(0,1,2,3)\) 的 DFT。\(n=4\)，\(\omega_4=i\)：

• \(y_0=0+1+2+3=6\)；
• \(y_1=0+1\cdot i+2\cdot i^2+3\cdot i^3=i-2-3i=-2-2i\)；
• \(y_2=0+1\cdot(-1)+2\cdot1+3\cdot(-1)=-2\)；
• \(y_3=0+1\cdot(-i)+2\cdot(-1)+3\cdot i=-2+2i\)。

所以 \(\mathrm{DFT}_4(0,1,2,3)=(6,\,-2-2i,\,-2,\,-2+2i)\)。

30.4.2 FFT：按奇偶下标分治

直接按 (30.8) 计算要 \(\Theta(n^2)\)。快速傅里叶变换（fast Fourier transform, FFT）利用单位复根的性质在 \(\Theta(n\lg n)\) 内算出 DFT。设 \(n\) 是 2 的幂（原书不讨论其他长度，见 30.7 节 chirp 变换）。

把系数按下标奇偶分成两组，定义两个次数界为 \(n/2\) 的多项式：

这样，"在 \(\omega_n^0,\dots,\omega_n^{n-1}\) 处求 \(A\)"化为：

1. 在 \((\omega_n^0)^2,(\omega_n^1)^2,\dots,(\omega_n^{n-1})^2\) 处求两个次数界 \(n/2\) 的多项式 \(A^{[0]}\)、\(A^{[1]}\)；
2. 按 (30.9) 合并。

由折半引理，这 \(n\) 个平方值只是 \(n/2\) 个 \(n/2\) 次单位根（各出现两次）。所以第 1 步恰好是两个规模减半、结构相同的 DFT 子问题。

推导拆解：用 \(n=4\)、\(a=(0,1,2,3)\) 走一遍，结果应与上面手算的 \((6,\,-2-2i,\,-2,\,-2+2i)\) 一致。

第一步，按奇偶拆分。\(A(x)=0+1x+2x^2+3x^3=(0+2x^2)+x(1+3x^2)\)，所以 \(A^{[0]}(y)=0+2y\)，\(A^{[1]}(y)=1+3y\)。这就是式 (30.9)：把 \(x^2\) 记作 \(y\)。

第二步，求值点平方。四个求值点 \(1,i,-1,-i\) 的平方是 \(1,-1,1,-1\)，所以 \(A^{[0]}\)、\(A^{[1]}\) 只需在 \(y=1\) 和 \(y=-1\) 两点求值（这是两个长度 2 的 DFT）：\(A^{[0]}(1)=2\)，\(A^{[0]}(-1)=-2\)，\(A^{[1]}(1)=4\)，\(A^{[1]}(-1)=-2\)。

第三步，按 (30.9) 合并。\(y_0=A(1)=A^{[0]}(1)+1\cdot A^{[1]}(1)=2+4=6\)；\(y_2=A(-1)=A^{[0]}(1)-1\cdot A^{[1]}(1)=2-4=-2\)；\(y_1=A(i)=A^{[0]}(-1)+i\cdot A^{[1]}(-1)=-2-2i\)；\(y_3=A(-i)=A^{[0]}(-1)-i\cdot A^{[1]}(-1)=-2+2i\)。

注意 \(y_0\) 和 \(y_2\) 共用同一对子结果，只差一个正负号；\(y_1\) 和 \(y_3\) 也是。这一对"一加一减"就是 30.6.1 节的蝴蝶操作。每一层都只做 \(n\) 次这样的加减，共 \(\lg n\) 层，所以总量是 \(n\lg n\)；而直接按 (30.8) 计算，每个 \(y_k\) 都要 \(n\) 次乘法，共 \(n^2\) 次。

import cmath

def recursive_fft(a):
    """CLRS RECURSIVE-FFT，约定 ω_n = e^{+2πi/n}，len(a) 为 2 的幂"""
    n = len(a)
    if n == 1:
        return [complex(a[0])]        # 单元素的 DFT 就是它自己：y0 = a0·ω_1^0
    wn = cmath.exp(2j * cmath.pi / n)
    w = 1
    y0 = recursive_fft(a[0::2])       # 偶下标系数
    y1 = recursive_fft(a[1::2])       # 奇下标系数
    y = [0j] * n
    for k in range(n // 2):
        t = w * y1[k]                 # 旋转因子 ω_n^k 乘 y1[k]，只算一次
        y[k] = y0[k] + t
        y[k + n // 2] = y0[k] - t
        w *= wn                       # 维护运行变量，不重复计算 ω_n^k
    return y

正确性：递归返回后，

因子 \(\omega_n^k\) 以正、负两种形式出现，称为旋转因子（twiddle factors）。

运行时间：除递归调用外，每层的 for 循环做 \(\Theta(n)\) 工作，

30.5 逆 DFT 与卷积定理

30.5.1 在单位复根处插值

求值的另一半是插值：已知 \(y=\mathrm{DFT}_n(a)\)，求 \(a\)。把 (30.8) 写成矩阵乘积 \(y=V_na\)：

定理 30.7 \(V_n^{-1}\) 的第 \((j,k)\) 元为 \(\omega_n^{-kj}/n\)。

证明：验证 \(V_n^{-1}V_n=I_n\)。其 \((j,j')\) 元为

白话解释：这个证明的思路是"猜一个答案，再验证它乘原矩阵等于单位矩阵"。之所以能猜中，是因为 \(V_n\) 的各列彼此"正交"：第 \(j\) 列与第 \(j'\) 列（取共轭后）对应元素相乘再相加，结果在 \(j=j'\) 时为 \(n\)，否则为 0。这与组合理论里"两个不相关因子的协方差为零"是同一种结构。各列正交的矩阵求逆很便宜：转置并取共轭，再除以 \(n\) 就行，不需要高斯消元。

以 \(n=4\)、\(j'-j=1\) 为例，求和引理说 \(1+i+i^2+i^3=1+i-1-i=0\)；\(j'-j=2\) 时是 \(1+(-1)+1+(-1)=0\)。这就是"等距分布在圆上的向量相加互相抵消"。（\(\overline{V_n}\) 上面的横线表示对每个元素取复共轭，即把 \(a+bi\) 换成 \(a-bi\)；"酉矩阵"是正交矩阵在复数里的对应物，乘它不改变向量长度，所以不放大误差。）

于是

把 (30.11) 与 (30.8) 对比：只要把 \(a\) 与 \(y\) 的角色互换、\(\omega_n\) 换成 \(\omega_n^{-1}\)、最后每项除以 \(n\)，就从 DFT 得到逆 DFT。所以逆 DFT 也能用 FFT 在 \(\Theta(n\lg n)\) 内算出（习题 30.2-4）。换句话说，\(V_n^{-1}=\frac1n\overline{V_n}\)：DFT 矩阵除以 \(\sqrt n\) 后是酉矩阵。这正是单位复根作插值点的数值优势——酉变换不放大误差，与 30.2.2 节提到的一般插值的病态性形成对照。

30.5.2 卷积定理

把求值和插值合起来，就得到多项式乘法的完整算法：

定理 30.8（卷积定理，convolution theorem） 对任意两个长度为 \(n\)（\(n\) 为 2 的幂）的向量 \(a,b\)，

30.5.3 补零与循环卷积：最常见的错误

为什么一定要补零到 \(2n\)？如果不补零，直接把两个长度 \(n\) 的 DFT 逐点相乘再逆变换，得到的是循环卷积（circular convolution）：

推导拆解：沿用前面的 \(a=(1,2)\)、\(b=(3,1)\)，正确的卷积是 \((3,7,2)\)。若不补零、直接用长度 2 的 DFT，只有两个求值点 \(\omega_2^0=1\) 和 \(\omega_2^1=-1\)。在这两个点上 \(x^2=1=x^0\)，所以 \(2x^2\) 和常数 2 无法区分。逆变换得到的是 \((3+2,\,7)=(5,7)\)：最高次项 2 被"绕回"加到了第 0 项上。按循环卷积公式核对：\(\tilde c_0=a_0b_0+a_1b_{(0-1)\bmod2}=3+2\cdot1=5\)，\(\tilde c_1=a_0b_1+a_1b_0=7\)。

补零到长度 4，即 \(a=(1,2,0,0)\)、\(b=(3,1,0,0)\)，乘积最高只到 \(x^2\)，四个求值点足以区分 \(x^0,x^1,x^2,x^3\)，就得到正确的 \((3,7,2,0)\)。

在时间序列里，循环卷积意味着"序列的开头混进了序列末尾的数据"——对回测来说，就是把未来数据混进了过去。30.8.1 节会用代码展示这一点。一般原则：长度 \(T\) 与 \(L\) 的两个序列做线性卷积，FFT 长度至少取 \(T+L-1\)。

30.6 高效 FFT 实现

30.6.1 蝴蝶操作

信号处理对速度要求极高，所以原书专门讨论如何把 FFT 写得更快。迭代版和递归版一样是 \(\Theta(n\lg n)\)，但常数可能更小（视实现而定，递归版有时对缓存更友好）。

第一个改进已经出现在上面的代码里：\(\omega_n^ky_k^{[1]}\) 在循环中要用两次（一加一减），称为公共子表达式（common subexpression），存到临时变量 \(t\) 里只算一次。"乘以旋转因子、存入 \(t\)、再与 \(y_k^{[0]}\) 分别相加和相减"这一组操作称为蝴蝶操作（butterfly operation，原书图 30.3），因为画成数据流图形如蝴蝶。

30.6.2 位逆序置换

观察 \(n=8\) 时的递归调用树（原书图 30.4）：

• 根：\((a_0,a_1,\dots,a_7)\)；
• 第二层：\((a_0,a_2,a_4,a_6)\) 与 \((a_1,a_3,a_5,a_7)\)；
• 第三层：\((a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)\)；
• 叶子顺序：\(a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7\)。

如果一开始就把输入排成叶子顺序，就可以自底向上计算：先两两做一次蝴蝶得到 \(n/2\) 个 2 元 DFT，再两两合并成 \(n/4\) 个 4 元 DFT……直到一个 \(n\) 元 DFT。

叶子顺序有一个漂亮的规律，称为位逆序置换（bit-reversal permutation）：\(a_k\) 应放在位置 \(\mathrm{rev}(k)\)，\(\mathrm{rev}(k)\) 是把 \(k\) 的 \(\lg n\) 位二进制表示倒过来读。0,4,2,6,1,5,3,7 的二进制是 000,100,010,110,001,101,011,111，倒过来恰好是 0,1,…,7。原因是：顶层按最低位（奇偶）分左右子树，下一层按次低位分，每层剥掉一位，最低位决定了最高层的走向。

30.6.3 迭代 FFT

import cmath

def bit_reverse_copy(a):
    n = len(a); L = n.bit_length() - 1
    A = [0j] * n
    for k in range(n):
        A[int(format(k, f'0{L}b')[::-1], 2)] = complex(a[k])
    return A

def iterative_fft(a, inverse=False):
    """CLRS ITERATIVE-FFT；inverse=True 时用 ω^{-1} 并除以 n，即逆 DFT"""
    n = len(a)
    A = bit_reverse_copy(a)
    sign = -1 if inverse else 1
    m = 2
    while m <= n:                         # 第 s 层，m = 2^s：合并成 m 元 DFT
        wm = cmath.exp(sign * 2j * cmath.pi / m)
        for k in range(0, n, m):          # 每组
            w = 1
            for j in range(m // 2):       # 组内 m/2 个蝴蝶
                t = w * A[k + j + m // 2]
                u = A[k + j]
                A[k + j] = u + t
                A[k + j + m // 2] = u - t
                w *= wm
        m *= 2
    return [x / n for x in A] if inverse else A

运行时间：BIT-REVERSE-COPY 每次反转 \(O(\lg n)\) 位，共 \(O(n\lg n)\)；实践中 \(n\) 事先已知，可查表做到 \(\Theta(n)\)，或用思考题 17-1 的摊还"逆序二进制计数器"。最内层循环的总执行次数

30.6.4 并行 FFT 电路

原书图 30.5 把迭代 FFT 画成电路：先做位逆序置换，然后是 \(\lg n\) 级（stage），每级 \(n/2\) 个蝴蝶并行执行。第 \(s\) 级有 \(n/2^s\) 组、每组 \(2^{s-1}\) 个蝴蝶，使用旋转因子 \(\omega_m^0,\dots,\omega_m^{m/2-1}\)（\(m=2^s\)）。电路深度（从输入到输出经过的计算元件最多有几层）为 \(\Theta(\lg n)\)，总共 \(\Theta(n\lg n)\) 个蝴蝶。这就是 GPU、FPGA 上 FFT 能做得极快的原因：每一级内部完全并行。

30.6.5 自己实现一遍，并与 numpy 对照

下面的脚本用上面两段代码验证 30.4 节的手算例子、位逆序顺序，并复现 30.1 节的多项式乘法。

import numpy as np
# recursive_fft, bit_reverse_copy, iterative_fft 定义同上

y = recursive_fft([0, 1, 2, 3])
print("DFT(0,1,2,3) =", np.round(y, 10))
print("numpy n*ifft  =", np.round(4 * np.fft.ifft([0, 1, 2, 3]), 10))
print("numpy fft     =", np.round(np.fft.fft([0, 1, 2, 3]), 10), "(符号约定相反，结果为共轭)")
print("bit-reverse order of 0..7:", [int(x.real) for x in bit_reverse_copy(list(range(8)))])

# A = 6x^3+7x^2-10x+9, B = -2x^3+4x-5（升幂系数）
a = [9, -10, 7, 6]; b = [-5, 4, 0, -2]
n = 8                                     # 2n，补零避免循环卷积混叠
ya = iterative_fft(a + [0] * (n - len(a)))
yb = iterative_fft(b + [0] * (n - len(b)))
c = iterative_fft([p * q for p, q in zip(ya, yb)], inverse=True)
print("A*B 系数(升幂):", np.round(np.real(c)).astype(int))
print("直接卷积      :", np.convolve(a, b))

rng = np.random.default_rng(0)
x = rng.standard_normal(1024)
err = np.max(np.abs(np.array(iterative_fft(list(x))) - 1024 * np.fft.ifft(x)))
print("n=1024 迭代版与 numpy 最大误差: %.2e" % err)

输出：

DFT(0,1,2,3) = [ 6.+0.j -2.-2.j -2.+0.j -2.+2.j]
numpy n*ifft  = [ 6.+0.j -2.-2.j -2.+0.j -2.+2.j]
numpy fft     = [ 6.+0.j -2.+2.j -2.+0.j -2.-2.j] (符号约定相反，结果为共轭)
bit-reverse order of 0..7: [0, 4, 2, 6, 1, 5, 3, 7]
A*B 系数(升幂): [-45  86 -75 -20  44 -14 -12   0]
直接卷积      : [-45  86 -75 -20  44 -14 -12]
n=1024 迭代版与 numpy 最大误差: 1.02e-12

两点要记住：

• 符号约定。原书 DFT 用 \(\omega_n=e^{+2\pi i/n}\)，numpy.fft.fft 用 \(e^{-2\pi i/n}\)。所以原书的 \(\mathrm{DFT}_n(a)\) 等于 n * np.fft.ifft(a)，而 np.fft.fft(a) 是它的复共轭（实输入时）。做卷积时两种约定都对，因为正逆变换成对使用；但在期权定价等需要具体积分方向的场合要对准公式。
• 工程实现。实际中用 numpy.fft/scipy.fft（底层为 pocketfft）或 FFTW。原书章末注记介绍：FFTW（"fastest Fourier transform in the West"）先运行一个"规划器"（planner）试跑几种分解方式，挑出在当前机器上最快的；它适配缓存，小规模子问题用展开的直线代码，并且对任意长度（包括大素数）都是 \(\Theta(n\lg n)\)。scipy.fft.next_fast_len 可以帮你挑一个因子只含 2、3、5 的"快长度"来补零。FFT 通常归功于 Cooley 与 Tukey（1960 年代），但 Heideman 等人考证出高斯在 1805 年就用过同样的思想。

30.7 思考题与习题选讲

原书本章的思考题大多是 FFT 思想的推广，其中几题在量化中直接有用。

集合的笛卡尔和计数（习题 30.1-7）：\(A\)、\(B\) 各含 \(n\) 个 \([0,10n]\) 内的整数，求所有 \(x+y\)（\(x\in A,y\in B\)）以及每个和出现的次数。把集合写成多项式 \(\sum_{a\in A}x^a\)，两者相乘，乘积中 \(x^s\) 的系数就是和为 \(s\) 的出现次数。用 FFT 在 \(O(n\lg n)\) 内完成。这是"用多项式乘法做计数"的标准套路，30.8.3 节的损失分布计算就是它的概率版本。

Karatsuba 分治乘法（思考题 30-1）：\((ax+b)(cx+d)\) 只需 3 次乘法：\(ac\)、\(bd\)、\((a+b)(c+d)\)，中间项为 \((a+b)(c+d)-ac-bd\)。递归下去得到 \(\Theta(n^{\lg3})\approx\Theta(n^{1.585})\) 的多项式乘法和大整数乘法。它比 FFT 慢，但没有浮点误差，中等规模时常数小。

Toeplitz 矩阵（思考题 30-2）：沿每条对角线为常数的矩阵（\(a_{ij}=a_{i-1,j-1}\)）称为 Toeplitz 矩阵。它只需 \(2n-1\) 个数表示；Toeplitz 矩阵乘向量可化为卷积，用 FFT 在 \(O(n\lg n)\) 内完成。量化意义：平稳时间序列 \((x_1,\dots,x_n)\) 的协方差矩阵 \(\Sigma_{ij}=\gamma(|i-j|)\) 正是对称 Toeplitz 矩阵，长序列的 GLS、高斯似然计算可借此加速（scipy.linalg.matmul_toeplitz、solve_toeplitz）。

多维 FFT（思考题 30-3）：\(d\) 维 DFT 可以沿每一维依次做一维 DFT，维的顺序无关，总时间 \(O(n\lg n)\)（\(n=n_1\cdots n_d\)），与维数无关。图像处理大量使用二维 FFT。

在一点求所有阶导数（思考题 30-4） 与 多点求值（思考题 30-5）：前者把 \(A(x_0+\omega_n^k)\) 的计算写成卷积，\(O(n\lg n)\) 得到全部 \(A^{(t)}(x_0)\)；后者构造"余数树"\(Q_{ij}=A\bmod\prod_{k=i}^j(x-x_k)\)，在 \(O(n\lg^2n)\) 内对任意 \(n\) 个点求值（习题 30.2-7 用同样的乘积树在 \(O(n\lg^2n)\) 内由根构造多项式）。

模算术 FFT / 数论变换（思考题 30-6）：复数 FFT 有舍入误差。若系数是整数，可以在模素数 \(p=kn+1\) 的整数环里做 FFT：取 \(\mathbb Z_p^*\) 的生成元 \(g\)，令 \(w=g^k\bmod p\) 作主 \(n\) 次单位根，DFT 与逆 DFT 在模 \(p\) 下互逆，结果精确。习题 30.2-6 是同一思想的另一个版本（模 \(2^{tn/2}+1\)、以 \(2^t\) 为单位根）。这需要第 31 章的数论知识。

chirp 变换 / Bluestein 算法（习题 30.2-8）：对任意复数 \(z\) 计算 \(y_k=\sum_ja_jz^{kj}\)。利用恒等式 \(kj=\frac{k^2+j^2-(k-j)^2}{2}\)，

30.8 量化实战

下面五段代码都只依赖 numpy/scipy/pandas，数据用随机模拟生成。

30.8.1 任意核的滚动加权统计

设收益序列长 \(T\)，要计算 \(s_t=\sum_{j=0}^{L-1}w_jr_{t-j}\)。权重 \(w\) 可以是截断的指数衰减、线性衰减、某个因子衰减曲线拟合出来的核，或者一个 FIR 滤波器。pandas 的 rolling().apply 对任意核是 \(\Theta(TL)\) 的 Python 循环，ewm 只支持指数核。FFT 卷积对任意核都是 \(\Theta((T+L)\lg(T+L))\)。

import numpy as np, time, pandas as pd
from scipy.signal import fftconvolve

rng = np.random.default_rng(42)
T = 200_000
r = 0.01 * rng.standard_normal(T)            # 模拟收益

L = 2000
lam = 0.999
w = (1 - lam) * lam ** np.arange(L)         # 截断的指数权重；可换成任意核
w /= w.sum()

def rolling_direct(r, w):
    L = len(w); out = np.full(len(r), np.nan)
    for t in range(L - 1, len(r)):
        out[t] = np.dot(w, r[t - L + 1:t + 1][::-1])
    return out

t0 = time.perf_counter(); s_dir = rolling_direct(r, w); t_dir = time.perf_counter() - t0
t0 = time.perf_counter()
s_fft = fftconvolve(r, w, mode="full")[:T]  # full 卷积的前 T 项就是因果滚动和
s_fft[:L - 1] = np.nan                       # 窗口未满的部分丢弃
t_fft = time.perf_counter() - t0
print(f"直接法 {t_dir:.2f}s, FFT 法 {t_fft:.3f}s, 加速 {t_dir / t_fft:.0f} 倍, "
      f"最大差 {np.nanmax(np.abs(s_dir - s_fft)):.1e}")

# 不补零的后果：循环卷积把序列尾部"绕回"开头
n = len(r)
circ = np.real(np.fft.ifft(np.fft.fft(r) * np.fft.fft(w, n)))   # 长度 n，没有补零
lin = fftconvolve(r, w, mode="full")[:n]
bad = np.abs(circ - lin)
print(f"未补零：前 L-1 个点最大误差 {bad[:L-1].max():.2e}；之后最大误差 {bad[L-1:].max():.1e}")
print(f"  第 0 个点：循环卷积 {circ[0]:+.2e}，线性卷积 {lin[0]:+.2e}（循环版混入了序列末尾 L-1 个未来收益）")

# 滚动方差 = 滚动均值(r^2) - 滚动均值(r)^2，两次卷积搞定
win = np.ones(250) / 250
m1 = fftconvolve(r, win)[:T]; m2 = fftconvolve(r ** 2, win)[:T]
var_fft = (m2 - m1 ** 2)[249:]
var_pd = pd.Series(r).rolling(250).var(ddof=0).to_numpy()[249:]
print(f"250 日滚动方差：FFT 与 pandas 最大相对差 {np.max(np.abs(var_fft / var_pd - 1)):.1e}")

输出：

直接法 0.30s, FFT 法 0.009s, 加速 34 倍, 最大差 8.7e-19
未补零：前 L-1 个点最大误差 1.56e-04；之后最大误差 1.1e-18
  第 0 个点：循环卷积 -1.42e-04，线性卷积 +3.52e-06（循环版混入了序列末尾 L-1 个未来收益）
250 日滚动方差：FFT 与 pandas 最大相对差 3.6e-14

解读：

• fftconvolve(r, w, mode="full") 内部自动补零到 \(T+L-1\) 以上，所以它的前 \(T\) 项就是只用过去数据的因果滚动和，没有前视。直接法即使用了向量化的 np.dot，也慢了几十倍（计时随机器和运行状态波动，重复运行在 30–80 倍之间）；若 \(L\) 更长、或要对几千只股票各算一遍，差距更大。
• 不补零时，前 \(L-1\) 个点被序列末尾的数据污染（误差量级与信号本身相当），之后的点恰好正确。这种 bug 很隐蔽：大部分数值对，只有开头一段错，而开头这一段在回测里恰恰用到了"未来"。
• "均值的平方差"公式在均值远大于标准差时会有严重的相消误差（例如直接对价格而不是收益算滚动方差）。对价格类序列，先减去一个参考值再算，或者改用 Welford 型的增量算法。

30.8.2 自相关函数与周期图

序列 \(x\) 的样本自相关在滞后 \(k\) 处是 \(\sum_tx_tx_{t+k}\)（归一化后）。它是 \(x\) 与自身翻转的卷积，所以也能用 FFT 计算：\(|\mathrm{DFT}(x)|^2\) 的逆变换就是（循环）自相关——这是 Wiener–Khinchin 关系的离散形式。同样要补零到至少 \(2N-1\) 才能得到线性自相关。\(|\mathrm{DFT}(x)|^2/N\) 本身叫周期图（periodogram），它在某个频率处的大小衡量序列在该频率上的"能量"，可用来检测周期性。ACF 的统计含义见第 06 册第 02a 章。

白话解释：DFT 可以理解为一次"回归分解"：把序列拆成各种频率的正弦、余弦波的叠加，\(|\mathrm{DFT}(x)|^2\) 在第 \(k\) 个频率上的值，大致相当于"周期为 \(N/k\) 的那条正弦波能解释多少方差"。所有频率上的周期图加起来，等于序列的总平方和（帕塞瓦尔恒等式），所以周期图就是把方差按频率做的分解。日内成交量有明显的"一天一个 U 形"，所以周期 48 根 K 线的那条波解释了大部分方差，在周期图上表现为一个尖峰。

自相关和周期图互为傅里叶变换（Wiener–Khinchin），所以两者携带的信息相同：一个从"滞后多少期还相关"的角度描述，一个从"哪个周期上能量多"的角度描述。\(|z|\) 表示复数 \(z=a+bi\) 的模长 \(\sqrt{a^2+b^2}\)。

下面模拟 60 个交易日、每日 48 根 5 分钟 K 线的成交量，带日内 U 形（开盘收盘放量）。

import numpy as np

rng = np.random.default_rng(7)
days, bars = 60, 48
tod = np.arange(bars)
u_shape = 1.0 + 1.5 * ((tod - bars / 2) / (bars / 2)) ** 2      # 开盘、收盘放量
vol = np.concatenate([u_shape * rng.lognormal(0, 0.35, bars) for _ in range(days)])
x = np.log(vol); x = x - x.mean()
N = len(x)

def acf_direct(x, K):
    d = np.dot(x, x)
    return np.array([np.dot(x[:N - k], x[k:]) / d for k in range(K + 1)])

def acf_fft(x, K):
    nfft = 1 << (2 * len(x) - 1).bit_length()          # 补零到 >= 2N-1，避免循环相关
    f = np.fft.rfft(x, nfft)
    ac = np.fft.irfft(f * np.conj(f), nfft)[:K + 1]
    return ac / ac[0]

K = 100
a1, a2 = acf_direct(x, K), acf_fft(x, K)
print("ACF 两种算法最大差: %.1e" % np.max(np.abs(a1 - a2)))
print("滞后 1, 24, 48, 96 的自相关:", np.round(a2[[1, 24, 48, 96]], 3))

P = np.abs(np.fft.rfft(x)) ** 2 / N                     # 周期图
freq = np.fft.rfftfreq(N, d=1.0)                         # 单位：每根 K 线的周期数
top = np.argsort(P[1:])[::-1][:3] + 1
for i in top:
    print(f"频率 {freq[i]:.5f} 周期 {1 / freq[i]:6.1f} 根K线  功率 {P[i]:.1f}")

输出：

ACF 两种算法最大差: 1.7e-16
滞后 1, 24, 48, 96 的自相关: [ 0.429 -0.419  0.415  0.411]
频率 0.02083 周期   48.0 根K线  功率 124.9
频率 0.04167 周期   24.0 根K线  功率 1.8
频率 0.50000 周期    2.0 根K线  功率 1.2

周期图在周期 48 根 K 线（一个交易日）处有一个比其他频率高两个数量级的峰；ACF 在滞后 48、96 处约 0.41，在滞后 24（半天，开盘对午盘）处为负。实务中，这种分析用来做成交量曲线预测（VWAP 执行的日内成交量分布，见第 07 册）和日内季节性调整：先用周期分量去掉日内模式，再研究残差。

注意：周期图在频率分辨率上受样本长度限制（相邻频点间隔 \(1/N\)），并且对噪声的方差不随 \(N\) 下降；正式的谱估计要做平滑（Welch 方法、多窗口法，scipy.signal.welch）。金融收益序列的周期性通常很弱，强峰多出现在成交量、波动率、价差这类有明显日内模式的变量上。

30.8.3 组合损失分布：多项式乘法的概率版本

设有 \(n\) 笔独立的信用敞口，第 \(m\) 笔以概率 \(p_m\) 违约、违约损失为整数单位 \(l_m\)。第 \(m\) 笔损失的概率生成函数是多项式 \(G_m(x)=(1-p_m)+p_mx^{l_m}\)；独立随机变量之和的生成函数是各自的乘积：

金融直觉：概率生成函数（probability generating function）是一个记账工具：把"损失为 \(s\) 的概率"写在 \(x^s\) 前面当系数。两笔贷款，损失分别为 1 和 2 单位、违约概率 0.1 和 0.2：

\[\big(0.9+0.1x\big)\big(0.8+0.2x^2\big)=0.72+0.08x+0.18x^2+0.02x^3.\]

逐项读：总损失 0 的概率 \(0.9\times0.8=0.72\)（都不违约）；损失 1 是 0.08（只有第一笔违约）；损失 2 是 0.18（只有第二笔违约）；损失 3 是 0.02（都违约）。多项式乘法自动把"指数相加"（损失相加）和"系数相乘"（独立事件概率相乘）配在一起，这就是独立随机变量之和的分布等于卷积的原因。

500 笔贷款就是 500 个这样的因子连乘。逐个相乘的方法 A 相当于每次做一个小卷积；方法 B 则利用"点值表示下乘法只是逐点相乘"：每个因子在单位复根处的值能直接写出，乘完后只做一次逆变换。拿到完整的损失分布后，VaR、ES（预期亏空）、任意分位数都可以直接读出，不需要蒙特卡洛模拟，也没有抽样误差。

import numpy as np

rng = np.random.default_rng(1)
n_names = 500
pd_ = rng.uniform(0.005, 0.05, n_names)          # 违约概率
lgd_units = rng.integers(1, 21, n_names)         # 违约损失，以 1 万元为单位，取整
M = int(lgd_units.sum()) + 1                     # 总损失可能取 0..M-1

# 方法 A：逐笔做两点卷积的递推，总 O(n*M)
dist = np.zeros(M); dist[0] = 1.0
for p, l in zip(pd_, lgd_units):
    new = dist * (1 - p)
    new[l:] += dist[:-l] * p
    dist = new

# 方法 B：每个因子在单位根处求值、相乘、一次逆 FFT
nfft = 1 << (M - 1).bit_length()                 # nfft >= M，不会绕回
k = np.arange(nfft)
logF = np.zeros(nfft, dtype=complex)
for p, l in zip(pd_, lgd_units):
    logF += np.log((1 - p) + p * np.exp(-2j * np.pi * k * l / nfft))
dist_fft = np.real(np.fft.ifft(np.exp(logF)))[:M]
print("两种方法最大差: %.1e" % np.max(np.abs(dist - dist_fft)))
cdf = np.cumsum(dist)
el = np.dot(np.arange(M), dist)
var99 = np.searchsorted(cdf, 0.99)
print(f"期望损失 {el:.1f} 万元（理论 {np.dot(pd_, lgd_units):.1f}），99% VaR = {var99} 万元")

输出：

两种方法最大差: 1.2e-16
期望损失 145.2 万元（理论 145.2），99% VaR = 258 万元

要点：

• 这里用的是 numpy 的负号约定，正变换的核 \(e^{-2\pi ikl/N}\) 与 np.fft.ifft 配对。
• 用对数相加再取指数，避免 500 个模长小于 1 的复数直接连乘时下溢；复对数的分支选择不影响 \(\exp(\sum\log)\) 的结果。
• FFT 长度必须不小于总损失的可能取值个数 \(M\)，否则高损失的概率会绕回到低损失上。
• 这个"独立违约"模型过于简单（真实组合有违约相关性，通常用单因子模型：先对系统因子条件化，条件独立下用本方法，再对因子积分），但卷积这一步的计算结构不变。保险精算中的聚合损失模型（Panjer 递推的替代）、两只骰子点数和之类的离散分布计算，都是同一个技巧。

30.8.4 频域滤波与前视偏差

"把价格做 FFT、去掉高频分量、再逆变换"是一种常见的去噪手法。问题是：全样本 FFT 在每个时刻的输出都用到了整段数据，包括未来。

import numpy as np
from scipy.signal import fftconvolve

rng = np.random.default_rng(1)
T = 1000
price = np.cumsum(rng.standard_normal(T + 50))
def lowpass_full_sample(x, keep=20):
    F = np.fft.rfft(x - x.mean()); F[keep:] = 0
    return np.fft.irfft(F, len(x)) + x.mean()
f1 = lowpass_full_sample(price[:T])          # 只用到 T 为止的数据
f2 = lowpass_full_sample(price[:T + 50])     # 多看了 50 个未来点
print(f"t=500 处滤波值：用前 {T} 点 {f1[500]:.3f}；加入 50 个未来点后 {f2[500]:.3f}")
print(f"t={T-1}（样本末端）处：{f1[T-1]:.3f} vs {f2[T-1]:.3f}，真实价格 {price[T-1]:.3f}")
# 因果替代：单边卷积核（只用过去数据），如截断 EWMA
w = 0.9 ** np.arange(60); w /= w.sum()
causal = fftconvolve(price[:T + 50], w)[:T + 50]
causal_short = fftconvolve(price[:T], w)[:T]
print("因果滤波：加入未来数据后前 T 点是否完全不变:", np.allclose(causal[:T], causal_short))

输出：

t=500 处滤波值：用前 1000 点 -18.597；加入 50 个未来点后 -17.903
t=999（样本末端）处：-26.342 vs -51.777，真实价格 -54.253
因果滤波：加入未来数据后前 T 点是否完全不变: True

两个现象：

1. 历史某一点（\(t=500\)）的滤波值会随着后面数据的到来而改变。在回测中用全样本滤波后的序列生成信号，等于在 \(t\) 时刻偷看了未来。
2. 样本末端的值严重失真（滤波值 \(-26.3\)，真实价格 \(-54.3\)）。原因是 DFT 默认序列是周期的，首尾不相接时会在边界产生大的伪影（Gibbs 现象）——而样本末端恰恰是实盘中唯一需要的那一点。

因果滤波器（单边卷积核）则没有这两个问题：加入未来数据后，过去的输出一个都不变。结论：FFT 可以放心地用来加速因果卷积，但不要把"全样本频域操作"的结果当作当时可得的信号。

30.8.5 FFT 期权定价（Carr–Madan）

若已知对数价格 \(\ln S_T\) 的特征函数 \(\phi(u)=E[e^{iu\ln S_T}]\)（Black–Scholes、Heston、Variance Gamma 等模型都有解析式），Carr 与 Madan（1999）给出用一次 FFT 计算整条执行价网格上欧式看涨期权价格的方法。

设对数执行价 \(k=\ln K\)，看涨价格为 \(C(k)\)。\(C(k)\) 在 \(k\to-\infty\) 时趋于 \(S_0\)，不可积；乘上阻尼因子 \(e^{\alpha k}\)（\(\alpha>0\)）后的 \(c(k)=e^{\alpha k}C(k)\) 可积，其傅里叶变换有解析式

白话解释：这一段的思路可以拆成三句话。

第一，特征函数是分布的另一种"身份证"。\(\phi(u)=E[e^{iu\ln S_T}]\) 和矩母函数 \(E[e^{t\ln S_T}]\) 形式相同，只是把实数 \(t\) 换成了虚数 \(iu\)；换成虚数的好处是它对任何分布都存在。正态分布 \(N(\mu,\sigma^2)\) 的特征函数是 \(e^{iu\mu-\sigma^2u^2/2}\)，代码中的 phi_bs 就是这个式子。很多模型（Heston、跳跃模型）写不出密度函数，却能写出特征函数，所以"从特征函数直接算期权价格"很有用。

第二，阻尼因子是为了让积分收敛。执行价趋于 0（\(k\to-\infty\)）时，看涨期权价格趋于 \(S_0\)，不会衰减到 0，对 \(k\) 积分会发散。乘上 \(e^{\alpha k}\) 后，左端被压到 0，就可以做傅里叶变换；算完再乘 \(e^{-\alpha k}\) 还原。

第三，为什么一次能算一整条期权链。\(C(k)\) 的公式是对 \(v\) 的积分，离散化后对每个执行价 \(k_u\) 都是一个和式 \(\sum_j(\cdots)e^{-ivjk_u}\)。选择网格使 \(\eta\lambda=2\pi/N\)，这 \(N\) 个和式就正好构成一个 DFT，用一次 FFT 全部算出。这和 30.4 节"一次算出多项式在 \(n\) 个点上的值"是同一件事。

import numpy as np
from scipy.stats import norm

S0, r, q, sigma, T = 100.0, 0.03, 0.0, 0.25, 0.5

def phi_bs(u):
    """Black–Scholes 下 ln S_T 的特征函数"""
    mu = np.log(S0) + (r - q - 0.5 * sigma ** 2) * T
    return np.exp(1j * u * mu - 0.5 * sigma ** 2 * u ** 2 * T)

def carr_madan(phi, N=4096, eta=0.25, alpha=1.5):
    lam = 2 * np.pi / (N * eta)                  # 对数执行价网格间距：η·λ = 2π/N
    b = N * lam / 2
    j = np.arange(N)
    v = eta * j
    k = -b + lam * j                             # 对数执行价网格
    psi = np.exp(-r * T) * phi(v - (alpha + 1) * 1j) / (alpha ** 2 + alpha - v ** 2 + 1j * (2 * alpha + 1) * v)
    w = eta / 3 * (3 + (-1) ** (j + 1)); w[0] = eta / 3   # Simpson 权重
    x = np.exp(1j * b * v) * psi * w
    C = np.exp(-alpha * k) / np.pi * np.real(np.fft.fft(x))  # 一次 FFT 得到整条执行价网格
    return np.exp(k), C

def bs_call(K):
    d1 = (np.log(S0 / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S0 * np.exp(-q * T) * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

K_grid, C_grid = carr_madan(phi_bs)
for K in [80, 90, 100, 110, 120]:
    c_fft = np.interp(np.log(K), np.log(K_grid), C_grid)
    print(f"K={K:3d}  FFT {c_fft:8.4f}   BS 解析 {bs_call(K):8.4f}")
mask = (K_grid > 50) & (K_grid < 200)
print("执行价 50~200 内网格点数 %d，最大绝对误差 %.1e" % (mask.sum(), np.max(np.abs(C_grid[mask] - bs_call(K_grid[mask])))))

输出：

K= 80  FFT  21.8351   BS 解析  21.8351
K= 90  FFT  13.7913   BS 解析  13.7908
K=100  FFT   7.7611   BS 解析   7.7603
K=110  FFT   3.8987   BS 解析   3.8986
K=120  FFT   1.7674   BS 解析   1.7669
执行价 50~200 内网格点数 226，最大绝对误差 2.2e-07

在 FFT 网格点上，价格与解析解相差不到 \(10^{-6}\)；表中个别执行价有 \(10^{-3}\) 量级的差，来自网格点之间的线性插值。这暴露了标准 FFT 的一个限制：\(\eta\lambda=2\pi/N\) 把积分网格和执行价网格绑在一起——想让执行价网格更密（\(\lambda\) 小），积分步长 \(\eta\) 就得变大，积分精度下降。解决办法正是 30.7 节的 chirp 变换（分数 FFT），它允许 \(\eta\lambda\) 任意选取。COS 方法、Lewis 公式属于同一思路。模型校准时要在几千组参数下反复计算整条期权链，FFT 方法的价值在此。期权与特征函数的背景见第 08 册。

本章小结

多项式有系数和点值两种表示：系数表示求值、加法快，乘法要 \(\Theta(n^2)\)；点值表示加法、乘法都是 \(\Theta(n)\)。在单位复根处求值就是 DFT，插值就是逆 DFT；单位复根的消去、折半、求和三条引理使得二者都能用分治的 FFT 在 \(\Theta(n\lg n)\) 内完成，从而卷积（多项式乘法）只需 \(\Theta(n\lg n)\)。实现上，迭代版 FFT 用位逆序置换加原地蝴蝶操作，并行电路深度为 \(\Theta(\lg n)\)。使用卷积定理时必须补零到足够长度，否则得到循环卷积——在时间序列中就是把未来数据绕回到开头。在量化中，FFT 用来加速任意核的滚动加权统计、自相关与谱分析、离散损失分布的卷积和基于特征函数的期权定价；但全样本的频域滤波会引入前视偏差，回测中只能用因果滤波。

练习

基础

1. 用式 (30.1)(30.2) 直接计算 \((7x^3-x^2+x-10)(8x^3-6x+3)\)，再按 30.2.4 节的四步方案用 FFT 重做一遍，核对结果。（原书 30.1-1、30.2-3。答案：\(56x^6-8x^5-34x^4-53x^3-9x^2+63x-30\)。）
2. 证明推论 30.4，并手算 \(\mathrm{DFT}_8\) 中 \(\omega_8\) 的各次幂（写成 \(a+bi\) 形式）。（原书 30.2-1。）
3. 写出逆 DFT 的伪代码，说明它与 FFT 的三处差别。（原书 30.2-4。提示：\(a,y\) 互换，\(\omega_n\to\omega_n^{-1}\)，除以 \(n\)。）
4. 跟踪 ITERATIVE-FFT 计算 \((0,2,3,-1,4,5,7,9)\) 的 DFT：写出位逆序后的数组，以及每一层结束后的数组。（原书 30.3-1。用 30.6.5 节的代码核对。）
5. 设集合 \(A=\{1,3,4\}\)、\(B=\{0,2,5\}\)。用多项式乘法求所有 \(a+b\) 及其出现次数。（原书 30.1-7 的小规模版本。）
6. 在 30.8.1 节的代码中，若把 fftconvolve(r, w, mode="full")[:T] 换成 mode="same"，结果会出什么问题？（提示：same 模式输出以核的中心对齐，相当于用了未来 \(L/2\) 期数据。）

进阶

7. 说明 Horner 法则如何改写成综合除法，在 \(\Theta(n)\) 时间内求 \(A(x)\) 除以 \((x-x_0)\) 的商 \(q(x)\) 与余数 \(r\)，并证明 \(r=A(x_0)\)。（原书 30.1-2。）
8. 设 \(n\) 是 3 的幂，推广 FFT 为三路分治，写出递归式并求解。（原书 30.2-5。提示：\(A(x)=A^{[0]}(x^3)+xA^{[1]}(x^3)+x^2A^{[2]}(x^3)\)，\(T(n)=3T(n/3)+\Theta(n)\)。）
9. 证明 Toeplitz 矩阵乘向量可在 \(O(n\lg n)\) 内完成，并用 scipy.linalg.matmul_toeplitz 验证：对 AR(1) 过程的协方差矩阵 \(\Sigma_{ij}=\phi^{|i-j|}/(1-\phi^2)\)，比较直接矩阵乘法与 Toeplitz 方法在 \(n=5000\) 时的速度。（原书思考题 30-2(c)。）
10. 用 chirp 变换实现任意长度 \(n\) 的 DFT，并在 \(n=1009\)（素数）上与 np.fft.fft 对比。（原书 30.2-8。提示：注意 numpy 的符号约定，\(z=e^{-2\pi i/n}\)；卷积长度取 \(\ge 2n-1\) 的 2 的幂。）
11. 在 30.8.3 节中引入单因子违约相关：设 \(Z\sim N(0,1)\)，条件违约概率 \(p_m(Z)=\Phi\big((\Phi^{-1}(p_m)-\sqrt\rho Z)/\sqrt{1-\rho}\big)\)。对 \(Z\) 用 Gauss–Hermite 积分，每个节点上用 FFT 求条件损失分布，再加权平均。比较 \(\rho=0\) 与 \(\rho=0.2\) 时的 99% VaR。
12. 用数论变换（思考题 30-6）在模 \(p=17\) 下计算 \((0,5,3,7,7,2,1,6)\) 的 DFT（取 \(g=3\)）。（提示：\(n=8\)，\(k=2\)，\(w=3^2\bmod17=9\)。）

原书推荐习题：30.2-2，30.1-7，30.2-8（chirp 变换），思考题 30-1（Karatsuba），思考题 30-2(c)（Toeplitz），思考题 30-6（数论变换）。

原书对照