本文框架性内容来自LiiiSTEM社区每周一题中相关内容 (特别鸣谢: 陈子豪). 本文以更自然的思路加以叙述, 同时对各步骤均作出详细证明, 考察相关操作的各项特性, 并在北太天元平台仿真验证, 提出算子改进方法.

🚀 Liii STEM - 世界最快数学输入法, 永久免费使用！用我的专属链接注册, 立享7天高阶OCR不限量免费使用: 👉 https://liiistem.cn/register.html?inviteCode=SEMLIMTK 或使用邀请码: SEMLIMTK, 期待你的加入.

摘要

傅里叶基适合描述循环系统, 有着信号学, 离散数学, 计算方法等不同领域的诸多原因.

本文通过递进的叙述路线,

从物理问题出发, 以环形阵列的邻值平均为切入点, 揭示邻值算数平均和循环移位操作之间的关系;

在线性代数层面, 系统分析变换结构, 通过证明邻值平均可分解为恒不变与循环移位的组合, 将问题转化为求循环移位算子的特征向量和特征值, 并发现循环移位的自然特征向量不是常见的坐标轴方向, 而是复指数振荡模式;

在傅里叶分析框架下, 经酉对角化从数字域变到频域, 将原变换解耦为傅里叶基上的独立缩放, 并重申了酉矩阵的优秀及酉性与能量守恒、数值稳定性和坐标独立性之间的关系.

背景

实验室里有一圈首尾相连的传感器, 均匀地排在一条环形轨道上. 每个传感器都会记录当前位置的某个信号强度, 比如温度, 噪声或者振动. 某一天, 老师拿到一组长度为 $N$ 的测量数据:

$$x = (x_0, x_1, \ldots, x_{N - 1})^T . $$

可是这组数据有些波动. 有的点突然偏高, 有的点突然偏低. 我们并不想立刻丢掉这些数据, 因为这些起伏里可能既有真实信息, 也有测量误差. 于是他提出了一个最朴素的处理办法: 做平均. 平均之后序列的第 $n$ 个数据定义为

$$y_n = \frac{x_n + x_{n + 1}}{2} . $$

因为传感器排成一圈, 所以最后一个点的右边并不是空着的, 而是又回到了第一个点. 因此更准确的表达应该写成

$$ y_n = \frac{x_n + x_{(n+1)\ \mathrm{mod}\ N}}{2},\quad n = 0,1,\ldots,N-1. $$

当 $n = N - 1$ 时,

$$y_{N - 1} = \frac{x_{N - 1} + x_0}{2} . $$

邻值平均变换的相关性质

线性

拿到一个变换首先需要考察其线性与否.

定义映射 $T : \mathbb{C}^N \rightarrow \mathbb{C}^N$ 为

$$ (T x)_n = \frac{x_n + x_{(n+1)\ \mathrm{mod}\ N}}{2},\quad n = 0,1,\ldots,N-1., n = 0, 1, \ldots, N - 1. $$

对于不同变元, 对任意常数 $c$, 对每个下标 $n$, 有

$$ (T (x + x^{\prime}))_n = \frac{(x_n + x_n^{\prime}) + (x_{n + 1} + x_{n + 1}^{\prime})}{2} = \frac{x_n + x_{n + 1}}{2} + \frac{x_n^{\prime} + x_{n + 1}^{\prime}}{2} = (T x)_n + (T x^{\prime})_n . $$

因此

$$T (x + x^{\prime}) = T x + T x^{\prime} . $$

同理,

$$ (T (c x))_n = \frac{c x_n + c x_{n + 1}}{2} = c \frac{x_n + x_{n + 1}}{2} = c (T x)_n . $$

所以

$$T (c x) = c T x. $$

故 $T$ 是线性变换.

矩阵表示

对于和向量, 序列有关的变换, 使用矩阵表示有利于研究.

邻点平均线性变换 $T$, 满足

$$ (T x)_n = \frac{x_n + x_{n + 1 \text{ {}} \text{ {}} (\operatorname{mod} N)}}{2} . $$

从公式本身看, $T$ 对第 $n$ 个位置做了三件事:

取当前位置自己的数值 $x_n$;
取右边邻居的数值 $x_{n + 1 \text{ {}} \text{ {}} (\operatorname{mod} N)}$;
将这两项取算术平均.

也就是说, 邻点平均可以看作两路信号的均值: 一条是原来的信号 $x$, 另一条是把整条信号整体向左平移后得到的信号.

因此定义一个线性变换 $S$ 以表示循环移位:

$$ (S x)_n = x_{n + 1 \text{ {}} \text{ {}} (\operatorname{mod} N)}, n = 0, 1, \ldots, N - 1. $$

例如, 当 $N = 4$ 时, 若

$$x = (x_0, x_1, x_2, x_3)^T, $$

则

$$S x = (x_1, x_2, x_3, x_0)^T . $$

于是, 邻点平均就可以写成

$$T x = \frac{1}{2} (x + S x), $$

可提取

$$T = \frac{1}{2} (I + S) . $$

在信号元素层面, 算数平均内容 = 输入向量和输入向量循环位移后结果的算术平均;

在矩阵 (线性变换) 层面, 算数平均变换 = 恒等变换与循环移位变换的算数平均.

寻找特征向量

矩阵作用于向量时, 可能在多个方向上同时拉伸, 旋转, 剪切, 各自由度之间彼此纠缠. 找到一组特殊的坐标轴 (即特征向量构成的基) 使得在这组新坐标下, 矩阵的作用退化为最简单的独立缩放 (每个特征向量方向只被乘上自己的特征值), 不同方向之间不再相互干扰.

对于一个线性变换, 想要寻找这组坐标就需要找到非零向量 $v$, 其在变换作用之后, 不会改变自己的基本形状, 而只是整体乘上一个常数. 也就是说, 我们希望找到 $v \neq 0$ 和某个复数 $\lambda$, 满足特征方程

$$S v = \lambda v. $$

如果某个信号满足这个关系, 那么 $S$ 在这个方向上的作用就被解耦了. 它不会再把该模式和别的模式混在一起, 而只是给它乘上一个常数.

在某个线性变换下只产生独立缩放的向量是最适合描述该线性变换的向量, 找到这个 (这些) 向量就能把矩阵进行基变换, 使得该线性变换变得单纯. 寻找特征向量的组合就是寻找基矩阵.

所以, 自然要问: 有没有某些信号, 在整体向左挪一格 (循环移位) 之后, 除了整体乘常数以外, 形状完全不变?

从直觉上来看:

常数信号 $(1, 1, 1, \ldots, 1)^T$ 被循环移位后完全不变;
信号 $(1, - 1, 1, - 1, \ldots \text{ })^T$ 被移位后会变成它的相反数, 但它的振荡形状并没有变;

因此我们猜想, 最适合描述移位的信号, 应该是那种每向前走一格, 就乘上同一个常数比例的模式.

于是考虑形如

$$v = (1, r, r^2, \ldots, r^{N - 1})^T $$

的等比数列形式向量. 由定义可知

$$S v = (r, r^2, \ldots, r^{N - 1}, 1)^T . $$

若满足

$$S v = \lambda v, $$

那么必须有

$$ (r, r^2, \ldots, r^{N - 1}, 1) = \lambda (1, r, r^2, \ldots, r^{N - 1}) . $$

比较第一个和最后一个分量, 得

$$ \left\{\begin{array}{l} \lambda = r\\ 1 = r \cdot r^{N - 1} = r^N \end{array}\right. . $$

因此, $r$ 必须满足

$$r^N = 1. $$

也就是说, $r$ 必须是 $N$ 次单位根.

联想到单位圆上不断旋转, 步长一致的向量, 设定步长 $\omega = e^{2 k \pi i / N}$, 其中 $k \in \mathbb{Z}, k \geqslant 0$, 则所有 $N$ 次单位根是

$$1, \omega^k, \omega^{2 k}, \ldots, \omega^{(N - 1) k} . $$

该序列集合可以充分描述所有的情况:

$N$ 表示可以回到原位的最少旋转次数
$k$ 表示旋转回到零相位时旋转过的圈数

注意到, $N$ 是和旋转结果有关的量, 而 $k$ 是和旋转过程有关的量. 一旋转向量旋转到第 $N$ 次 (或第某次) 时, 可能已旋转不止一圈; 因为我们仅关注旋转结果 (静态的位置), 所以这种周期性带来了信息的冗余.

进一步讲, $N$ 是将 $k$ 圈 ($2 k \pi$ 作为总相位) 进行划分. 若 $k \geqslant N$, 则每一步旋转过的相位 $2 k \pi / N$ 超过 $2 \pi$. 显然可以将该步实际旋转的相位对 $2 \pi$ 取模来得到其等效 1 圈内旋转的程度; 而该步实际旋转的相位 $2 k \pi / N$ 对 $2 \pi$ 取模即为 $k \operatorname{mod} (N)$, 因此可以通过取模或限制 $k$ 的范围 (令 $k = 0, 1, \ldots, N - 1$) 来避免剩余周期的无效讨论.

于是我们得到 $N$ 个特征向量:

$$ v_k = (1, \omega^k, \omega^{2 k}, \ldots, \omega^{(N - 1) k})^T, k = 0, 1, \ldots, N - 1. $$

该向量组各元素满足

$$S v_k = \omega^k v_k . $$

这说明对循环移位来说, 最自然的模式并不是通常的坐标轴方向, 而是这些按复指数振荡的模式.

特征值分解和相似对角化

多项式方程组的变量和系数可以分别被提取成对应矩阵, 因此我们可以将方程组写成矩阵乘法形式.

从代数角度, 我们可以分别观察研究系数矩阵和变量矩阵, 将矩阵的相关性质和方程组的结果, 类型相匹配;
从几何角度, 我们可以将方程形象化为线性的变换, 考察其对于变元向量的图形化 (拉伸旋转) 作用;
从应用角度, 矩阵可以通过并行化提高计算速度.

因此我们希望将

$$ S v_k = \omega^k v_k $$

这一方程组也写成矩阵形式.

对于 $v_k$, 自然想到从 $k$ 和 $N$ 两个维度展开, 把它们依次排列, 写成矩阵形式, 即令

$$V = \left[ v_0 \text{~} v_1 \text{~} \cdots \text{~} v_{N - 1} \right], $$

对于我们希望写成的矩阵形式特征方程的左边, 矩阵左乘单个列向量在几何层面可表示线性变换作用于该向量的结果.

矩阵左乘列向量组也可以表示该矩阵独立作用于各个列向量后再将作用结果组成矩阵的结果. 因为 $S V$ 的第一列就是 $S v_0$, 与其他列无关; 同样, 结果的每一列也都等同于 $S$ 矩阵左乘对应 $V$ 矩阵的那一列.

需要说明的是, 尽管将 $V$ 写成分块行向量的形式左乘 $S$ (即 $S \left[ v_0 \text{~} v_1 \text{~} \cdots \text{~} v_{N - 1} \right]$) 看上去维度并不匹配, $S$ 和 $V$ 矩阵从实际的元素展开来看都是 $N \times N$ 的, 所以 $S V$ 运算合法.

对于我们希望写成的矩阵形式特征方程的右边, 每一列应为 $\omega^k v_k$, 整体应为 $\left[ \omega^0 v_0 \text{~} \omega^1 v_1 \text{~} \cdots \text{~} w^{N - 1} v_{N - 1} \right]$. 可以通过构造对角矩阵, 左 (右) 乘实现行 (列) 向量缩放, 对应初等列 (行) 变换. 该对角矩阵的对角元素应为每一列向量所需的缩放倍数 $w^k$, 即

$$ D = \operatorname{diag}(1, \omega, \omega^2, \ldots, \omega^{N - 1}) = \left[ \begin{array}{cccc} 1 & 0 & \cdots & 0\\ 0 & \omega & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \omega^{N - 1} \end{array} \right], $$

这样特征方程组就可以写成矩阵形式:

$$S V = V D. $$

因为 $\omega^k$ ($k = 0, \ldots, N - 1$) 是 $N$ 个互不相同的 $N$次单位根，它们对应的特征向量 $v_k$ 线性无关，所以 $V$ 是可逆方阵, 即可得到标准的相似对角化

$$S = V D V^{- 1}.$$

酉对角化

酉对角化过程

酉矩阵要求列向量 (或行向量) 构成一组标准正交基, 而目前已有的 $V$ 矩阵满足正交性, 还差一个标准性. 注意到 $V$的每一元素模长都为 1, 因此 $V$ 的每一列向量的模长为 $\sqrt{N}$, 自然想到给 $V$矩阵乘上 $\sqrt{N}$. 令

$$ F_N = \frac{1}{\sqrt{N}} D = \frac{1}{\sqrt{N}} \left[ v_0 \text{~} v_1 \text{~} \cdots \text{~} v_{N - 1} \right] = \frac{1}{\sqrt{N}} (\omega^{j k})_{, 0 \leq j, k \leq N - 1} . $$

此即离散傅立叶矩阵.

于是

$$ S = \left( \sqrt{N} F_N \right) D \left( \sqrt{N} F_N \right)^{- 1} = \sqrt{N} F_N \hspace{0.17em} D \hspace{0.17em} \frac{1}{\sqrt{N}} \hspace{0.17em} F_N^{- 1} = F_N D F_N^{- 1}, $$

$F_N$ 仍满足相似对角化基矩阵的要求.

不仅如此, 我们还观察到 $F_N$ 的第 $k$ 列是 $u_k = v_k / \sqrt{N}$. 计算第 $k$ 列与第 $l$ 列的内积:

$$ {(F^{\ast}_{ N}}^{} F_N)_{k, l} = \langle u_k, u_l \rangle = \frac{1}{N} \langle v_k, v_l \rangle = \frac{1}{N} \sum_{j = 0}^{N - 1} (\omega^{j k}) \overline{(\omega^{j k}) } = \frac{1}{N} \sum_{j = 0}^{N - 1} \omega^{j (k - l)} . $$

• 当 $k = l$ 时, 每一项都是 1 ，求和得 $(F_N^{\ast} F_N)_{k, k} = 1$.

• 当 $k \neq l$ 时, 这是公比为 $ \omega^{k - l},\omega^{k - l} \neq 1 $ 的等比数列, 故:

$$ \sum_{j = 0}^{N - 1} \omega^{j (k - l)} = \sum_{j = 0}^{N - 1} (\omega^{k - l})^j = \frac{1 - (\omega^{k - l})^N}{1 - \omega^{k - l}} = \frac{1 - (\omega^N)^{k - l}}{1 - \omega^{k - l}} = \frac{1 - 1^N}{1 - \omega^{k - l}} = 0, $$

因此

$$F_N^{\ast} F_N = I, $$

即 $F_N$ 是酉矩阵, $S$ 可从相似对角化分解"进化"成酉对角化分解:

$$S = V D V^{- 1} = F_N D F_N^{- 1} = F_N D F_N^{\ast}. $$

在傅立叶基下, 循环移位变换 $S$ 不再把各个位置混在一起, 而是变成对角矩阵; 每一个傅立叶模式 $v_k$ 单独演化, 只会乘上自己的特征值 $\omega^k$.

因为

$$T = \frac{1}{2} (I + S), $$

所以

$$T = F_N \left( \frac{1}{2} (I + D) \right) F_N^{\ast} . $$

记

$$ \left\{\begin{array}{l} \Lambda = \frac{1}{2} (I + D) = \operatorname{diag} (\lambda_0, \lambda_1, ..., \lambda_{N - 1})\\ \lambda_k = \frac{1 + \omega^k}{2}, k = 0, 1, ..., N - 1, \end{array}\right. $$

那么

$$T = F_N \Lambda F_N^{\ast} . $$

可以看到, 原来位置上 (数字域) 的算数平均操作 $T$ 在傅立叶坐标下变成了每个复指数振荡模式按照 $\Lambda$ 单独缩放. 这就是为什么傅立叶分解在信号与系统中如此经典.

为什么需要酉对角化

避免求逆

矩阵求逆是数值计算中最危险的操作之一, 执行高斯消元或类似的求逆算法会引入舍入误差, 且条件数可能放大误差; 而 $F_N^{\ast}$ 只是共轭转置, 几乎没有计算代价, 没有任何数值不稳定的风险.

标准性, 能量守恒, 保范数保内积

酉矩阵最好的性质是保持向量的长度和内积不变, 意味着信号时频域变换不改变总能量. 这也对应 Parseval’s theorem:

$$\sum_{j = 0}^{N - 1} | x_j |^2 = \sum_{k = 0}^{N - 1} | \hat{x}_k |^2 . $$

如果 $F_N$ 不是酉矩阵 (比如相似对角化时的 $V$ 矩阵), 而只是非标准化的可逆矩阵, 那么变换前后的能量比例就被扭曲, 频域系数就不再与原信号的能量直接对应 (虽然也没什么大事只是比例问题).

正交性, 坐标独立

信号 $x$ 在 $u_k$ 基上的线性展开为:

$$x = \hat{x}_0 u_0 + \hat{x}_1 u_1 + \ldots + \hat{x}_{N - 1} u_{N - 1} . $$

其中 $U = (u_0 u_1 \ldots u_{N - 1})$.

若 $u_k$ 为一般基, 则$\hat{x} = U^{- 1} x$, 此时每个坐标分量都依赖于所有基向量;
若 $u_k$ 为酉基, 则 $U^{- 1} = U^{\ast}$, 逆矩阵退化成共轭转置, $\hat{x} = U^{\ast} x$, $\hat{x}_k = \langle x, u_k \rangle$, 每个坐标分量退化成简单的内积投影.

傅里叶基下, $\hat{x} = F_N^{\ast} x$. 此即频谱系数向量. 酉矩阵的正交性保证频域的每个分量都是原信号在对应振荡模式上的独立投影.

对称性, 提高算法效率

由于 $F_N^{- 1} = F_N^{\ast}$，DFT 与 IDFT 共享计算结构, 只差一个共轭 (甚至 $\omega^k$ 是纯虚数) 和归一化因子; 这也是 FFT 算法能够同时高效计算傅里叶正反变换的原因。

回顾离散傅里叶矩阵的产生思路:

原问题是邻点平均;
邻点平均的核心动作是循环移位;
为了看懂循环移位, 要寻找它的特征向量;
循环移位的特征向量恰好是复指数振荡模式的旋转向量;
把这些特征向量排成矩阵并归一化, 就得到傅立叶矩阵.

多次平均的低通效果

变换域 (频域) 效果

将多次平均操作

$$ T^m x = (F_N \Lambda F_N^{})^m x = F_N \Lambda^m F_N^{} x = F_N \Lambda^m \hat{x} . $$

按列展开写成频域叠加形式:

$$T^m x = \sum_{k = 0}^{N - 1} \lambda_k^m \hat{x}_k u_k, $$

其中 $u_k$ 是 $F_N$ 的第 $k$ 列, 即归一化特征向量 $u_k = \frac{1}{\sqrt{N}} (1, \omega^k, \omega^{2 k}, \ldots, \omega^{(N - 1) k})^T$.

利用欧拉公式,

$$1 + e^{i \theta} = 2 \cos (\theta / 2) e^{i \theta / 2}, $$

令 $\theta_k = 2 \pi k / N$, 得

$$ \lambda_k = \frac{1 + e^{i \theta_k}}{2} = \cos \hspace{-0.17em} \left( \frac{\pi k}{N} \right) e^{i \pi k / N}, $$

所以不同特征值的模长

$$| \lambda_k | = \left| \cos \left( \frac{\pi k}{N} \right) \right| . $$

可见 $\lambda_k$ 对 $k = 0, 1, ..., N - 1$ 单调递减.

由于对所有 $k \neq 0$ 都有 $| \lambda_k | < 1$.

当 $m \to \infty$ 时,

$$\lambda_k^m \to 0, $$

故

$$ \Lambda^m = \operatorname{diag}(\lambda_0^m, \lambda_1^m, ..., \lambda_{N - 1}^m) \longrightarrow \operatorname{diag}(1, 0, ..., 0) . $$

所以从变换域来看, $T^m$ 相当于对信号反复低通滤波:

每一轮平均都以衰减因子 $| \cos (\pi k / N) |^m$ 削弱不同交流分量,
直流分量因为 $\lambda_0 = 1$ 被完全保留,
当迭代次数趋于无穷，所有空间差异被抹平, 系统收敛到平衡态, 对应常数场.

空间域 (数字域) 效果

将以上结果代回谱分解:

$$ \lim_{m \to \infty} T^m = F_N \Lambda^m F_N^{\ast} = u_0 {u^{\ast}_0}^{} . $$

注意到 $u_0 = \frac{1}{\sqrt{N}} (1, 1, \ldots, 1)^T$, 故

$$ u_0 {u^{\ast}_0}^{} = \frac{1}{N} \left[ \begin{array}{cccc} 1 & 1 & \cdots & 1\\ 1 & 1 & \cdots & 1\\ \vdots & \vdots & \ddots & \vdots\\ 1 & 1 & \cdots & 1 \end{array} \right] . $$

对任意输入信号 $x$, 极限为

$$ \lim_{m \to \infty} T^m x = \left( \frac{1}{N} \sum_{j = 0}^{N - 1} x_j \right) \left[ \begin{array}{c} 1\\ 1\\ \vdots\\ 1 \end{array} \right] . $$

迭代 $m$ 次后, 当前位置的结果实际上是周围 $m$ 个邻居的加权平均, 其权重系数遵循二项式分布, 在 $m$ 较大时逼近高斯分布. 所以 $T^m x$ 的极限等于信号 $x$ 的算术平均值均摊到所有位置.

高阶取平均的操作实际上是一个扩散过程. 在连续极限下, 该算子对应于一维热传导方程:

$$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} . $$

$t$ 随着 $m$ 的增加, 高频信号 (热点, 噪声) 被快速抹平, 系统趋于均匀; 可据此感性理解算子的低通特性.

仿真验证

相位偏移

代码构造完成后, 我发现使用算子 $T$ 滤波后, 信号变平滑, 但发生了明显的水平漂移. 这种相位偏移显而易见是由 $T^m$ 变换中的 $k, m, N$ 共同决定的.

算子 $T$ 的特征值 (频率响应) 为:

$$\lambda_k = \frac{1}{2} (1 + e^{i \frac{2 \pi k}{N}}) $$

考虑 $m$ 次迭代系统总频率响应:

$$ H (k) = \lambda_k^m = \underbrace{\cos^m \left( \frac{\pi k}{N} \right)}_{\text{幅值响应}} \cdot \underbrace{e^{i \frac{\pi km}{N}}}_{\text{相位响应}} $$

即对应的相位滞后量为 $\pi k / N$.

根据傅里叶变换的位移性质, 频域的线性相移

$$e^{i \frac{2 \pi k}{N} \cdot \Delta n} $$

对应数字域平移量

$$\Delta n = \frac{m}{2}, $$

即信号在数字域上向左平移了 $m / 2$ 个采样点.

在 $N = 1024, k = 2$ 的条件下, 信号自身的周期为 $T_{sig} = N / k = 512$ 点. 当 $m = 512$ 时, 位移量 $\Delta n = 256$, 相对于信号自身周期的偏移率为

$$\frac{\Delta n}{T_{sig}} = \frac{256}{512} = \frac{1}{2} . $$

偏移半个周期意味着波峰与波谷完全对调, 故视觉上呈现反相.

不过以上只考虑了整数点数的数字域平移. 对于$\Delta n \not\in \mathbb{Z}$ 没有原生意义; 然而借助 DFT 的周期延拓, 非整数平移可以存在, 通过

DFT 将信号嵌入到由 $e^{i 2 \pi kn / N}$ 张成的有限维空间中.

考虑纯相位滤波器

$$H (k) = e^{i \frac{2 \pi k}{N} \Delta n}, k = 0, 1, ..., N - 1. $$

根据 DFT 时移性质，该相位因子在频域的作用等价于时域平移 $\Delta n$.

对输入信号频谱 $\hat{x}_k$ 逐点乘以 $H (k)$ 后做 IDFT 即可得到分数延迟后的输出:

$$ y [n] = \frac{1}{N} \sum_{k = 0}^{N - 1} \left( \hat{x}_k \hspace{0.17em} e^{i \frac{2 \pi k}{N} \Delta n} \right) \hspace{0.17em} e^{i \frac{2 \pi k}{N} n} = \sum_{m = 0}^{N - 1} x [m] \hspace{0.17em} h [n - m], $$

其中 $h [n]$ 是 $H (k)$ 的 IDFT，即该滤波器的时域脉冲响应.

$$ h [n] = \frac{1}{N} \sum_{k = 0}^{N - 1} e^{i \frac{2 \pi k}{N} (n + \Delta n)} . $$

令 $\theta = \frac{2 \pi}{N} (n + \Delta n)$, 利用等比数列求和:

$$ h [n] = \frac{1}{N} \sum_{k = 0}^{N - 1} e^{i k \theta} = \frac{1}{N} \cdot \frac{1 - e^{\mathrm{i} N \theta}}{1 - e^{i \theta}} = \frac{1}{N} \cdot \frac{1 - e^{i 2 \pi (n + \Delta n)}}{1 - e^{i \frac{2 \pi}{N} (n + \Delta n)}} . $$

利用欧拉公式化简分子分母:

$$ h [n] \hspace{0.17em} = \hspace{0.17em} \frac{1}{N} \hspace{0.17em} e^{i \pi (n + \Delta n) \hspace{0.17em} \frac{N - 1}{N}} \cdot \frac{\sin \left[ \pi (n \hspace{0.17em} + \hspace{0.17em} \Delta n) \right] }{\sin \left[ \frac{\pi}{N} (n \hspace{0.17em} + \hspace{0.17em} \Delta n) \right] } . $$

最后将频率索引对称化, 把 $k > N / 2$ 的项映射为负频率 $k - N$, 相位因子 $e^{i \pi (n + \Delta n) \frac{N - 1}{N}}$ 就被吸收掉了, 得到实值形式 (在共轭对称条件下):

$$ h [n] = \frac{\sin [\pi (n + \Delta n)]}{N \sin [\frac{\pi}{N} (n + \Delta n)]} ~ . $$

此式即为 Dirichlet 核, 也称周期 sinc (很熟悉, 之前 Gibbs 现象推导收敛的时候用到了).

双边优化

为了消除单边算子带来的相位漂移, 可以引入双边 (对称) 算子, 同时参考左右邻居:

$$T_{sym} = \frac{1}{4} S^{- 1} + \frac{1}{2} I + \frac{1}{4} S $$

其在信号域的表现为: $y_n = \frac{1}{4} x_{n - 1} + \frac{1}{2} x_n + \frac{1}{4} x_{n + 1}$.

其特征值为:

$$ \lambda_{k, sym} = \frac{1}{4} e^{- i \frac{2 \pi k}{N}} + \frac{1}{2} + \frac{1}{4} e^{i \frac{2 \pi k}{N}}, $$

因为

$$\cos \theta = \frac{e^{i \theta} + e^{- i \theta}}{2}, $$

所以:

$$ \lambda_{k, sym} = \frac{1}{2} + \frac{1}{2} \cos \left( \frac{2 \pi k}{N} \right) = \cos^2 \left( \frac{\pi k}{N} \right) $$

由于特征值 $\lambda_{k, sym}$ 是纯实数, 其相位始终为 0; 所以无论迭代多少次 $m$, 信号只会在原位被平滑.

附录仿真代码

MATLAB

clear; clc;
% ==========================================
% 单边算子 vs 双边(对称)算子：时频域对比
% 注意：请确保已在北太天元中开启 fft 插件
% ==========================================
N = 1024;      
n_idx = (0:N-1)';
k = 16;

% 1. 构造信号：2个周期的低频主信号 + 高频噪声
clean_signal = sin(2 * pi * k * n_idx / N); 
noise = 0.5 * randn(N, 1);                  
x = clean_signal + noise;

% 2. 执行 FFT 进入频域 (保留原始频谱以供对比)
X_hat = fft(x);
k_idx = (0:N-1)';

% 3. 构造算子频域响应 (特征值)
lambda_asym = 0.5 * (1 + exp(1i * 2 * pi * k_idx / N)); % 单边
lambda_sym  = 0.5 * (1 + cos(2 * pi * k_idx / N));      % 双边(对称)

% 4. 计算四个场景的频域系数
Y_asym_512  = X_hat .* (lambda_asym .^ 512);
Y_asym_2048 = X_hat .* (lambda_asym .^ 2048);
Y_sym_512   = X_hat .* (lambda_sym .^ 512);
Y_sym_2048  = X_hat .* (lambda_sym .^ 2048);

% 5. IFFT 转回时域 (取实部消除浮点误差)
y_asym_512  = real(ifft(Y_asym_512));
y_asym_2048 = real(ifft(Y_asym_2048));
y_sym_512   = real(ifft(Y_sym_512));
y_sym_2048  = real(ifft(Y_sym_2048));

% ===================== 绘图可视化 =====================
% 创建一个更大的窗口以容纳大字号文本
figure('Name', '时频域全景对比', 'Position', [50, 50, 1100, 950]);

% 将数据打包进元胞数组
y_time = {y_asym_512, y_asym_2048, y_sym_512, y_sym_2048};
Y_freq = {Y_asym_512, Y_asym_2048, Y_sym_512, Y_sym_2048};
row_titles = {
    '1. 单边算子 m=512',
    '2. 单边算子 m=2048',
    '3. 双边算子 m=512',
    '4. 双边算子 m=2048'
};

% 提取原始频谱的幅值 (0 到 N/2)
k_half = 0:N/2;
X_mag = abs(X_hat(1:N/2+1));

% --- 统一的字号设置参数 ---
ax_font = 16;     % 坐标轴数字字号
lbl_font = 16;    % 坐标轴标签字号
tit_font = 18;    % 标题字号
leg_font = 14;    % 图例字号

for i = 1:4
    % --- 左列：时域对比 ---
    subplot(4, 2, 2*i - 1);
    plot(n_idx, clean_signal, 'Color', [0.7 0.7 0.7], 'LineWidth', 1.5, 'LineStyle', '--'); hold on;
    plot(n_idx, y_time{i}, 'r-', 'LineWidth', 1.5);
    
    set(gca, 'FontSize', ax_font); % 调整坐标轴数字字号
    title([row_titles{i}, ' - 时域'], 'FontSize', tit_font, 'FontWeight', 'bold');
    ylabel('振幅', 'FontSize', lbl_font);
    grid on; xlim([0, N-1]); ylim([-1.5, 1.5]);
    
    if i == 4 
        xlabel('采样点 n', 'FontSize', lbl_font); 
    end
    if i == 1 
        legend('纯净信号', '滤波结果', 'Location', 'northeast', 'FontSize', leg_font); 
    end

    % --- 右列：频域频谱对比 ---
    subplot(4, 2, 2*i);
    Y_mag = abs(Y_freq{i}(1:N/2+1));
    
    plot(k_half, X_mag, 'Color', [0.7 0.7 0.7], 'LineWidth', 1); hold on;
    plot(k_half, Y_mag, 'b-', 'LineWidth', 1.5);
    
    set(gca, 'FontSize', ax_font); % 调整坐标轴数字字号
    title([row_titles{i}, ' - 频域幅值谱'], 'FontSize', tit_font, 'FontWeight', 'bold');
    ylabel('幅值', 'FontSize', lbl_font);
    grid on; 
    xlim([0, 50]); 
    
    if i == 4
        xlabel('频率指数 k', 'FontSize', lbl_font); 
    end
    if i == 1 
        legend('原始带噪频谱', '滤波后频谱', 'Location', 'northeast', 'FontSize', leg_font); 
    end
end
点击展开查看更多

摘要

背景