第一次被 Toeplitz 矩阵卡住，是几年前跑信道估计仿真的时候。信道矩阵正好是 Toeplitz 结构，当时用最朴素的循环按 O(n²) 硬乘，矩阵一上 2048 阶就慢得没法看。后来才搞清楚，这类矩阵天生就是给 FFT 准备的——把 Toeplitz 矩阵塞进循环矩阵，用快速傅里叶变换把矩阵向量乘法从 O(n²) 降到 O(n log n)。这个思路在科研仿真、嵌入式实时处理、FPGA 加速里都很常见，而且一旦理解透了，矩阵乘法的行观点和列观点、卷积和滤波、以及各种 FFT 库的选型，都会串成一条线。

这篇文章会从 Toeplitz 矩阵最基本的对角线结构讲起，然后拆解怎么构造循环矩阵、为什么要取前 m 行、行观点列观点分别在什么场景下有用，最后给出可以直接跑的 NumPy 实现，以及 STM32F4、Vivado FFT IP 核落地时最容易踩的坑。

1. Toeplitz 矩阵的“内存节省”和“计算浪费”并存

1.1 对角线相同的定义谁都会背，关键是理解它省掉了什么

Toeplitz 矩阵的定义一句话就能说完：每条对角线上的元素全部相同。

一个 m×n 的 Toeplitz 矩阵，虽然看起来有 m×n 个元素，但真正独立的只有 m+n-1 个。它的第一列和第一行决定一切。举个例子：

T = [ r0 r1 r2 ] [ c1 r0 r1 ] [ c2 c1 r0 ]

这里的c = [c0, c1, c2]是第一列，r = [r0, r1, r2]是第一行，其中c0 = r0是共同的左上角元素。

换成人话：老式居民楼每一层的户型都一样，只要记住一层楼的图纸，整栋楼长什么样就确定了。Toeplitz 矩阵也一样，第一列加第一行就把整个矩阵钉死了。

正因为这样，存储一个 1000×1000 的 Toeplitz 矩阵只需要大约 1999 个元素，而不是 100 万个。内存上是巨大优势，但很多人忽略了另一面——计算时如果还按照普通稠密矩阵去做乘法，那这份结构优势就被浪费得一干二净。

1.2 卷积、信道辨识、位移不变系统全都绕不开它

Toeplitz 矩阵不是数学家的玩具，它出现在所有“位移不变”的场景里。

最典型的就是卷积。一个 FIR 滤波器对信号做卷积，写成矩阵形式必然出现 Toeplitz 结构。通道均衡、匹配滤波、波束形成、地震数据处理、光学成像去模糊，凡是本身带平移不变性的问题，矩阵形式都是 Toeplitz 或分块 Toeplitz。

我做信道估计的时候，发射训练序列，接收信号和高斯噪声的关系写下来，中间那个信道卷积矩阵就是 Toeplitz。你要解它、要拿它做矩阵向量乘法，绕不开。

还有一种来源是相关矩阵。信号的自相关矩阵在很多条件下也能写成 Toeplitz，尤其是平稳随机过程。Wiener 滤波、AR 模型参数估计、谱估计里全是这玩意儿。

所以这不是一个“偶尔碰到”的矩阵，而是一个覆盖面极广的矩阵结构。

1.3 O(n²) 不是不能跑，是反复调用时真的顶不住

如果只做一次矩阵向量乘法，n=4096，O(n²) 大约是 1600 万次乘加，单精度浮点下现代 CPU 也就几毫秒的事，看起来完全能忍。

问题是真实场景里几乎没人只乘一次。

比如迭代求解 Toeplitz 线性方程组，一次迭代就要乘一次；做系统辨识，每个时间窗口都要对训练矩阵来一次；做图像复原，每次迭代都要对被模糊核构造的 Toeplitz 块矩阵做乘法。一个算法跑几百上千次迭代，O(n²) 立刻变成瓶颈。

更不用说 n 到 10 万量级的情况。O(n²) 是 100 亿次乘加，就算 CPU 再快也顶不住；而 O(n log n) 大概是 170 万次左右，这是天壤之别。

所以核心问题变成了：能不能利用 Toeplitz 矩阵的重复结构，把每次乘法的复杂度拉低？

答案就是 FFT。但直接拿 FFT 去乘 Toeplitz 矩阵是不行的，中间需要一层转换。

2. 算法核心：把 Toeplitz 塞进循环矩阵再交给 FFT

2.1 循环矩阵和 DFT 是一对天生搭档

要说清楚 Toeplitz 怎么用 FFT 加速，必须先讲它的近亲——循环矩阵。

循环矩阵长这样：

C = [ k0 k1 k2 k3 ] [ k3 k0 k1 k2 ] [ k2 k3 k0 k1 ] [ k1 k2 k3 k0 ]

每一行都是上一行循环右移一位。循环矩阵最漂亮的性质是：它可以用离散傅里叶变换对角化。

准确地说，任何循环矩阵 C 都可以写成：

C = F^{-1} · diag(FFT(k)) · F

其中 k 是循环矩阵的第一列，F 是 DFT 矩阵，FFT(k) 是 k 的离散傅里叶变换。于是 C 乘以任意向量 x 就可以这样算：

C x = IFFT( FFT(k) ⊙ FFT(x) )

这里 ⊙ 表示逐元素相乘。

一次循环矩阵向量乘法，只需要做两次 FFT 加一次 IFFT，再加上一次逐元素乘法，总复杂度 O(L log L)，比 O(L²) 快一个量级以上。

问题是 Toeplitz 矩阵并不是循环矩阵，它右上角和左下角不对称。但我们可以把 Toeplitz 矩阵“嵌入”到一个更大的循环矩阵里，让循环矩阵的一部分行和 Toeplitz 矩阵完全重合。

2.2 第一列怎么拼？一张小表格说清楚

假设 Toeplitz 矩阵 T 是 m×n 的，第一列为 c（长度 m），第一行为 r（长度 n），且 c[0] = r[0]。

构造一个长度为 L = m + n - 1 的序列 k，作为大循环矩阵的第一列：

k = [ c0, c1, ..., c_{m-1}, r_{n-1}, r_{n-2}, ..., r_1 ]

注意这里 r 的顺序是从尾到头，而且去掉了 r[0]，因为 r[0] 和 c[0] 是同一个元素，已经放在头部了。

再用一个具体的 3×3 例子来看：

c = [a0, a1, a2] r = [a0, a_{-1}, a_{-2}] T = [ a0 a_{-1} a_{-2} ] [ a1 a0 a_{-1} ] [ a2 a1 a0 ]

按上面的规则拼出：

k = [a0, a1, a2, a_{-2}, a_{-1}]

然后把 x 向量补零到同样长度：

x_pad = [x0, x1, x2, 0, 0]

现在建立以 k 为第一列的循环矩阵 C，维度 L×L。C 的前 m 行和前 n 个有效列，恰好就是原来的 Toeplitz 矩阵。

这一步就是整个加速算法最关键的地方，也是最容易出错的地方。我见过很多人要么忘了把 r 反转，要么把 r[0] 重复拼接了一次，结果构造出来的根本不是原矩阵。

2.3 为什么结果取前 m 行就刚刚好

先别急着写代码。要理解循环嵌入后为什么结果直接取前 m 行。

还是用 3×3 的例子展开。x_pad = [x0, x1, x2, 0, 0]，循环矩阵 C 的第一列是 k = [a0, a1, a2, a_{-2}, a_{-1}]。根据循环矩阵定义，C 的第 s 行元素是 k[(s-j) mod L]。

逐行算：

• 第 0 行： [a0, a_{-1}, a_{-2}, a2, a1]，和 x_pad 做点积得到 a0·x0 + a_{-1}·x1 + a_{-2}·x2，这正好是 T 第 0 行乘以 x。
• 第 1 行： [a1, a0, a_{-1}, a_{-2}, a2]，点积得到 a1·x0 + a0·x1 + a_{-1}·x2，正好是 T 第 1 行乘以 x。
• 第 2 行： [a2, a1, a0, a_{-1}, a_{-2}]，点积得到 a2·x0 + a1·x1 + a0·x2，正好是 T 第 2 行乘以 x。

看到规律了吗？循环矩阵的行与行之间不断循环移位，而 Toeplitz 矩阵的每一行也在往同一个方向滑动。只要 k 的第一列拼接方向正确，在循环卷积的前 m 行内，每一行的循环移位和 Toeplitz 的行结构完美对齐。

所以完整的算法流程是：

1. 根据 Toeplitz 的第一列 c 和第一行 r，构造 k，长度 L。
2. 把 x 补零到长度 L。
3. 计算 IFFT(FFT(k) ⊙ FFT(x_pad))。
4. 结果取前 m 个元素，就是 T·x。

看起来像在做循环卷积，实际上是借用循环矩阵的行移位结构，把 Toeplitz 的乘法“伪装”了一次。

3. 行观点和列观点在 Toeplitz 加速中的实际分工

3.1 行观点：每一行都是同一个模板被切开再平移

矩阵乘法里有个经典的说法：结果矩阵的每一行，等于左矩阵那一行和右矩阵整个矩阵做线性组合；每一列，等于左矩阵整个矩阵和右矩阵那一列做线性组合。

在 Toeplitz 矩阵里，行观点的意义特别直白。

T 的每一行其实都是从同一个“模板”上切下来的片段。第一行是模板从左边开始的 n 个元素，第二行是模板向右错一位的 n 个元素。整个 T 就是模板在不同位置上的快照。

所以 T·x 的每一行，本质上就是在做“模板和 x 的滑窗内积”。这就是相关运算的标准形式。而相关运算又能通过翻转和卷积联系起来，卷积又能用 FFT 加速。这条链路就是 Toeplitz 加速的理论底气。

行观点的价值在于：它让你意识到 T·x 不是在算 m×n 个独立的点积，而是在算同一个模板和 x 在 m 个不同位置上的匹配。既然模板是同一个，那模板本身的信息就能在频域里被压缩成一组固定的频域系数，一次算好，后面反复用。

3.2 列观点：批量矩阵乘法是靠按列拆分来复用的

如果你要算的不是 T·x，而是 T·X，其中 X 是一个 n×k 的普通矩阵，那该怎么办？

列观点直接给出答案：把 X 拆成 k 列，每一列单独做一次 Toeplitz 向量乘法，然后把结果按列拼回去。

T · X = [ T·x1, T·x2, ..., T·xk ]

这看起来像废话，但它引出了一个极其重要的工程优化点：k 次乘法中，FFT(k) 只需要计算一次。

每次 Toeplitz 向量乘法本来要做三次 FFT 规模的运算，分别是 FFT(k)、FFT(x)、IFFT(结果)。但 FFT(k) 只依赖 Toeplitz 矩阵本身，跟 x 无关。批量场景下，先算一次 FFT(k)，然后每一列只需要做一次 FFT 和一次 IFFT。

如果 X 有 1000 列，那就是从 3000 次级 O(L log L) 运算降到 2001 次级，省了接近三分之一。更重要的是，很多数值库支持批量 FFT，把一整个二维数组按 axis=0 一次性变换，比循环调用 1000 次单列 FFT 快得多，因为内存访问更连续，循环展开和向量化也更容易。

3.3 遇到转置和伴随，别让存储顺序坑了你

Toeplitz 矩阵的转置仍然是 Toeplitz 矩阵，但第一列和第一行会交换位置。

如果要用快速算法算 Tᵀ·X，不能直接把原来的 c、r 拿过来用，要先交换：

new_c = [r0, r1, ..., r_{n-1}] new_r = [c0, c1, ..., c_{m-1}]

而且注意 new_r 去掉 r0 后的反向拼接顺序也要跟着变。

实际写代码时还有存储顺序的问题。NumPy 默认是行优先，如果你的 X 是按行存的，想按列拆开做批量 FFT，直接切片 X[:, j] 会得到不连续的内存视图，性能反而差。更好的做法是把 X 转成列优先数组（Fortran order），或者用np.asfortranarray(X)，这样按列访问时 cache 友好很多。

这一点在大矩阵乘法里差距非常明显，我在第 4.2 节会再强调一次。

4. 最小可复现实现：NumPy 二十行搞定

4.1 先写一版能跑通的核心函数

理论说再多，不如直接看代码。下面这个函数同时支持 T·x 和 T·X 两种形式。

import numpy as np def toeplitz_mult_mat(c, r, X): c = np.asarray(c, dtype=float) r = np.asarray(r, dtype=float) m = len(c) n = len(r) # 理论最小长度是 m + n - 1，这里取 nextpow2 方便 FFT L = 1 while L < m + n - 1: L <<= 1 # 拼出循环矩阵第一列：c 正序 + r 除 r0 外的反序 k = np.concatenate([c, r[-1:0:-1]]) X = np.atleast_2d(X) if X.shape[0] != n: if X.shape[1] == n: X = X.T else: raise ValueError("X 的行数必须等于 Toeplitz 矩阵的列数 n") # 预计算 FFT(k)，批量场景下只算这一次 kF = np.fft.fft(k, L) # X 补零到长度 L，同时按列批量 FFT Xp = np.zeros((L, X.shape[1]), dtype=complex) Xp[:n, :] = X # 频域点乘，再做批量 IFFT Y = np.fft.ifft(kF[:, None] * np.fft.fft(Xp, axis=0), axis=0) # 取前 m 行就是结果 return Y[:m, :].real

这段代码里几个关键点：

• r[-1:0:-1]表示从 r 的最后一个元素开始，反向取到 r[1]，不包含 r[0]，正好是 n-1 个元素。
• kF[:, None]把一维向量变成列向量，方便和二维批量 FFT 结果逐列相乘。
• np.fft.fft(Xp, axis=0)一次处理所有列，而不是 for 循环逐列算，性能差距非常大。

4.2 长度为什么是 nextpow2(m+n-1)

这个坑我一开始踩得很深。

理论上循环嵌入只需要 L = m + n - 1 就能完整表示 Toeplitz 矩阵的全部行。但如果直接用这个长度做 FFT，FFT 算法本身也支持任意长度，NumPy 会用 Bluestein 算法处理非 2 的幂，慢不说，还会引入额外误差。

取 nextpow2 的好处有两个：

1. 快速傅里叶变换在长度是 2 的幂时最省时间。
2. 补零到更大的长度相当于增加了循环周期，保证循环卷积的回绕项不会污染前 m 行的结果。

有人会问：补零到 m+n-1 不就够了吗，为什么还要 nextpow2？

理论上是够了，但用 nextpow2 完全不影响结果的正确性，因为多出来的那些位置本来就不参与前 m 行的计算。反而是少取长度会出大问题——如果 L 小于 m+n-1，循环矩阵第一列里根本放不下完整的 r 反转序列，结果必然错。

所以稳妥做法永远是先保证 L ≥ m+n-1，再取 nextpow2。

4.3 用一个随机小矩阵验证正确性

写数值代码最忌讳直接上大矩阵。先用小矩阵验证，确保索引和方向都没问题。

from scipy.linalg import toeplitz c = np.array([1.0, 2.0, 3.0, 4.0]) r = np.array([1.0, 0.5, 0.25, 0.125]) T = toeplitz(c, r) X = np.random.randn(4, 6) Y_fast = toeplitz_mult_mat(c, r, X) Y_ref = T @ X print(np.max(np.abs(Y_fast - Y_ref)))

我实际跑出来的结果是1.4e-14左右，双精度下完全可接受。如果你看到结果差到 1e-1 量级，不用怀疑，一定是k的拼接方向错了，或者补零位置错了。

验证的时候建议同时检查：

• 单个向量的场景X = np.random.randn(4)。
• 多列矩阵场景X = np.random.randn(4, 6)。
• 非方阵场景，比如 m=4, n=6，或者 m=6, n=4。

非方阵最考验索引逻辑，因为很多人写代码时脑子默认是方阵，一旦 m≠n 就翻车。

4.4 单精度和双精度的误差实测

我把同一个测试在不同浮点精度下跑了几组，结果很有意思：

单精度误差看起来不大，但要注意这个误差会随着变换长度变大而缓慢增长，而且如果 Toeplitz 矩阵本身条件数很大，误差还会放大很多。做频谱分析、滤波这类对精度不敏感的任务，float32 完全够用；但要做高精度数值计算，比如解 Toeplitz 线性方程组，还是老老实实用 float64。

如果是在 STM32F4 这类单精度平台上做，建议先用 PC 上的 float32 仿真完整流程，确认误差在自己能接受的范围内，再移植到嵌入式。

5. FFT 选型和硬件落地：从 STM32F4 到 Vivado IP 核

5.1 FFT 为什么能到 O(n log n)：蝶形分解一句话版本

DFT 的朴素定义是：

X[k] = Σ x[n] · e^{-i 2π k n / N}

直接按这个公式算，每个频点要 N 次复数乘法，总共 N² 次。FFT 的核心是 Cooley-Tukey 分解——把长度为 N 的序列按奇偶分成两半，每个长度 N/2 的 DFT 再继续拆。每拆一层，运算量减半；拆 log₂N 层，总运算是 O(N log N)。

蝴蝶运算是这个分解在代码层面的体现。两个输入，一个旋转因子，算出两个输出。这一层换一层，理解上不需要太深，够用就行。

工程上真正要注意的不是 FFT 原理，而是你手头平台上的 FFT 实现怎么选、怎么调参数。

5.2 STM32F4 上用 CMSIS-DSP 的注意事项

STM32F4 系列带 FPU 和 DSP 指令，跑 FFT 很常见。CMSIS-DSP 库里有现成的arm_rfft_fast_f32（实数 FFT）、arm_cfft_f32（复数 FFT）、arm_cfft_q15/arm_cfft_q31（定点 FFT）。

实际使用中这几个注意点非常关键：

第一，内存对齐。CMSIS-DSP 的 FFT 函数对输入数组有对齐要求，特别是arm_cfft_f32要求数组按 4 字节对齐。定义全局数组基本没问题，但如果是在函数内部动态分配或栈上定义局部数组，最好加上__ALIGNED(4)或者直接用静态数组。

第二，单精度误差。STM32F4 的 FPU 是单精度的，4096 点 FFT 的误差大概 1e-5 量级，做频谱分析没问题。如果你要做的是 Toeplitz 矩阵相关的数值计算，这个误差可能不够，需要评估。

第三，预计算卷积核。如果要在 STM32 上做 Toeplitz 矩阵向量乘法，FFT(k)这部分完全可以预先算好，存成常量数组放在 Flash 里。每次实时处理只需要FFT(x)、逐元素复数乘法、IFFT。这样每次耗时就从三次 FFT 变成两次，省掉一次。

第四，实数信号优化。如果输入信号是实数，别用复数 FFT。用arm_rfft_fast_f32比arm_cfft_f32快差不多一半，因为它利用实数序列频谱的共轭对称性，只需要算一半频点。

5.3 Vivado FFT IP 核该选哪种架构

FPGA 平台上的 FFT 通常直接用 Xilinx 的 FFT IP 核，Vivado 里配置界面很直观，但几个选项背后的权衡要清楚。

FFT IP 核主要架构有三种：

如果做 Toeplitz 矩阵相关的高速卷积，肯定是 Pipelined Streaming 最合适，因为数据可以源源不断喂进去，输出连续。缺点是非常吃 DSP slice 和 BRAM。如果 FPGA 资源紧张，用 Radix-4 Burst I/O，一次处理一块数据，处理完再接收下一块，中间有间隙，但资源开销小很多。

配置 IP 时要注意 scaling schedule。Xilinx FFT IP 有几种缩放模式：

• Unscaled：不做缩放，大信号中间过程容易溢出。
• Scaled：每级手动配置缩放因子，需要自己算。
• Block Floating Point：自动缩放，输出带指数，精度有保障。

我个人建议大多数场景直接用 Block Floating Point。它会在运算过程中自动检测溢出并做缩放，输出的指数带你记录缩放倍数，误差比 Unscaled 小得多，又不用手动算每级缩放。缺点是多了一个指数输出通道，多占一点逻辑。

另外 AXI4-Stream 接口的握手必须处理对。tvalid 和 tready 信号时序搞错，IP 核直接卡死，仿真看不出来，上板就出问题。建议先在 Vivado 自带的 example design 上跑通，再改自己的数据通路。

6. 我踩过的几个坑，希望你能直接绕开

6.1 拼接序列的方向反了：小矩阵手工展开才发现

我第一次把算法从论文落到代码时，循环矩阵第一列 k 用的是np.concatenate([c, r[1:]])，也就是 r 正序拼接。小矩阵一验证，只有第一行是对的，下面全错。

原因前面已经说了：r 必须反转拼接，而且要排除 r[0]，因为 r[0] 就是 c[0]。

当时我检查了很久没看出问题，后来拿 3×3 矩阵手写逐行展开才定位到。所以我现在写这类代码，第一步永远是先拿 4×4 的矩阵，用纸笔把 k 写出来，再运行代码对比。这一步千万别省。

6.2 变换长度取少了一个值：循环混叠会让后几行彻底崩掉

另一个高频错误是 L 取成了max(m, n)或者m + n - 2。

当 L 小于理论最小长度m + n - 1时，循环矩阵的第一列装不下完整的反转 r，Toeplitz 矩阵右下角那些元素的贡献会被“回绕”到别的位置，表现出来就是后几行严重错误，而且不是均匀的误差，是直接算错。

这个问题最难排查的地方在于：小矩阵因为后续补零多，有时碰巧误差不明显；一旦 m、n 稍微大一点，结果立刻全碎。

解决很简单：先算出理论最小值，再往上取 2 的幂。

6.3 直接取 .real 之前要想清楚复数是怎么来的

FFT 和 IFFT 的中间结果全是复数。虽然 Toeplitz 矩阵本身是实矩阵，输入 x 也是实数，FFT 乘完再 IFFT 之后，理论上虚部应该接近 0，可以直接.real取值。

但要注意：

• 浮点误差会让虚部出现 1e-15 量级的残留值，直接忽略没问题。
• 如果你不小心把某个频谱的共轭对称结构弄坏了，比如只用了前半段频点做 IFFT，结果虚部会非常大，这时直接.real会丢掉真实信息。

我遇到过有人优化性能，改用np.fft.rfft处理实数序列，但忘了 k 的 FFT 也要满足共轭对称，结果数据乱套。要优化可以，把两边的对称性都处理好，再切到rfft/irfft。

6.4 并行批量处理时共享预计算频谱的隐患

批量处理 T·X 时，多线程并行每个列向量的 FFT 是可以的，kF是只读的，所有线程共享没问题。

真正的坑是 Xp 补零数组。如果每个线程各自写自己的那一列，互不干扰，没问题；但如果有人为了省内存，让多个线程共享整个 Xp 数组，某个线程在补零阶段可能覆盖到其他线程的数据。

这种事在 OpenMP 里很好发生。我给出的建议是：要么每个线程维护自己的列切片，要么把批量 FFT 交给 NumPy 一次性完成，别手动分线程去操作同一个二维数组的列。

最后想分享的一点体会

做 Toeplitz 相关计算这几年，我最大的感受是：FFT 只是加速的“下半场”，“上半场”是把矩阵结构吃透。

很多资料一上来就丢 FFT 公式和复杂度结论，但真正写代码时，最容易错的却不是 FFT 本身，而是 Toeplitz 怎么嵌入循环矩阵、第一列怎么拼、结果取哪几行。这些索引问题靠死记硬背记不住，最好的办法是先画一个小矩阵，手工展开算一遍，再上代码。

我现在的习惯是：不管项目多着急，先写一个 4×4 或 8×8 的参考实现验证索引，再谈性能优化。这个习惯帮我省下的调试时间，远比我花在验证上的几分钟多得多。后面有机会，我还会单独写一篇 Toeplitz 线性方程组的迭代求解，涉及预处理和 Levinson-Durbin 递推，那又是另一个有意思的话题了。