痛苦系列 | DSP-02 频域切片:DTFT与DFT的探索
文中代码及相关仿真脚本均可在 BlogCode 中获取,本篇代码集中于 DSP/dsp02 目录下。
痛苦系列 DSP 专栏 :抛弃死记硬背与纯符号堆砌,从工程中的真实痛点切入,剖析直觉误区,由第一性原理推导数学工具与物理边界,并结合 Python 仿真验证。
前言
关于痛苦的延续
上一篇我们在连续与离散的交界处划了一刀,讨论了采样定理、时域原子序列与线性时不变系统(LTI)的卷积和。写完那篇后,我不禁回忆起求学时期学习医学影像处理时的一个疑惑:既然时域卷积已经能够描述系统的动态响应,为什么又出现了傅里叶变换?
当年在校学习时,专业课程的教材编排不够系统——既缺乏从实际物理问题出发的层层引导,也缺少严谨连贯的公式推导链条,许多变换与核心概念几乎是从天而降、突兀登场。这种“知其公式却不知其所以然”的困惑,伴随了我很长时间。
直到后来深入从事脑机接口(BCI)、经皮神经电刺激等神经工程实践,天天与真实的生理电信号打交道,我才深刻体会到频域工具的不可替代性:在实际工程中,我们采集到的传感器数据,在时域上往往就是一串充斥着工频干扰、肌电漂移和电刺激伪迹的离散数值序列。如果你只戴着“时域眼镜”,面对这团混杂的数据,想要用时域差分或者阈值条件把高频毛刺剥离出来,几乎等同于在面粉里挑出特定粒径的沙子——不仅破坏了原始波形的时序特征,更常常因为算法的局部震荡导致后续特征提取全面失真。
为了看清信号的内在骨架,我们必须转换观察坐标系。本篇我们将从数学推导与物理约束出发,拆解离散时间傅里叶变换(DTFT)与离散傅里叶变换(DFT)的底层逻辑,并结合十个硬核的 Python 仿真实验,彻底搞清楚频域分析的威力与边界。
阅读导图
为方便读者诸君把握脉络,本篇将聚焦于以下几个核心问题:
DTFT 的本质与物理约束 :从正交投影到实序列共轭对称性,频域解耦的数学证明与物理意义。
卷积定理推导与工程价值 :为什么时域滑动卷积能等价于频域标量相乘?
从连续走向离散的桥梁——DFT :频域采样到底带来了什么?时域周期延拓的本质。
能量守恒与坐标系变换 :帕斯瓦尔定理的严密验证,以及归一化因子 1 / N 1/N 1/ N 的物理来源。
工程幻象避坑指南 :栅栏效应与频谱泄漏的内在机理、窗函数的取舍代价、时域补零(Zero-padding)的插值真相。
最隐蔽的工程陷阱——循环卷积 :时域周期折叠混叠的数学根源,以及恢复线性卷积的补零边界条件。
物理极限与时宽带宽积 :Dirichlet 核的演化、时频不确定性的物理约束。
一、 DTFT:理想频域世界的正交投影与物理约束
在上一篇中,我们把连续信号切片为离散时间序列 x[n]。对于一个绝对可和的时域序列(即 ∑ n = − ∞ ∞ ∣ x [ n ] ∣ < ∞ \sum_{n=-\infty}^\infty |x[n]| < \infty ∑ n = − ∞ ∞ ∣ x [ n ] ∣ < ∞ ),我们定义其**离散时间傅里叶变换(DTFT)**为:
X ( e j ω ) = ∑ n = − ∞ ∞ x [ n ] e − j ω n X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n} X ( e j ω ) = n = − ∞ ∑ ∞ x [ n ] e − j ω n
这里的 ω \omega ω 被称为数字角频率 ,量纲是 rad/sample \text{rad/sample} rad/sample (弧度/样本)。
如果我们把模拟连续信号的采样周期记为 T s T_s T s ,采样频率记为 f s = 1 / T s f_s = 1/T_s f s = 1/ T s ,连续模拟角频率记为 Ω = 2 π f \Omega = 2\pi f Ω = 2 π f 。那么离散序列中的采样点 n n n 对应物理时间 t = n T s t = n T_s t = n T s 。将模拟复指数信号 e j Ω t e^{j\Omega t} e j Ω t 代入采样点,得到:
e j Ω n T s = e j ( Ω T s ) n = e j ω n e^{j\Omega n T_s} = e^{j (\Omega T_s) n} = e^{j\omega n} e j Ω n T s = e j ( Ω T s ) n = e j ω n
由此可知数字频率与模拟物理频率之间的换算关系:
ω = Ω T s = 2 π f f s \omega = \Omega T_s = 2\pi \frac{f}{f_s} ω = Ω T s = 2 π f s f
当模拟频率达到奈奎斯特极限 f = f s / 2 f = f_s / 2 f = f s /2 时,数字角频率恰好达到 ω = π \omega = \pi ω = π 。因此,数字角频率的区间 [ − π , π ] [-\pi, \pi] [ − π , π ] (或 [ 0 , 2 π ] [0, 2\pi] [ 0 , 2 π ] )完整映射了连续信号在采样后不发生混叠的全部物理频段。
从几何投影的角度来看,式中的 e − j ω n e^{-j\omega n} e − j ω n 是一组连续的正交复指数基底。DTFT 本质上就是在计算时域序列 x [ n ] x[n] x [ n ] 在不同频率基底上的复内积。内积模值的大小,代表了该频率分量在信号中所占的权重。
1. 单频正弦信号的 DTFT 映射
我们来看第一个仿真实验:一个纯净的单频离散正弦信号 x [ n ] = sin ( ω 0 n ) x[n] = \sin(\omega_0 n) x [ n ] = sin ( ω 0 n ) 。
在实验参数中,设物理采样率 f s = 100 Hz f_s = 100\text{ Hz} f s = 100 Hz ,信号物理频率 f 0 = 10 Hz f_0 = 10\text{ Hz} f 0 = 10 Hz 。此时数字角频率为:
ω 0 = 2 π 10 100 = 0.2 π ≈ 0.628 rad/sample \omega_0 = 2\pi \frac{10}{100} = 0.2\pi \approx 0.628\text{ rad/sample} ω 0 = 2 π 100 10 = 0.2 π ≈ 0.628 rad/sample
(图1:单频正弦信号的时域波形、DTFT 幅度谱、相位谱及复平面轨迹)
在图 1 中,左上角是离散时间波形 x [ n ] x[n] x [ n ] ;右上角是其对应的 DTFT 幅度谱 ∣ X ( e j ω ) ∣ |X(e^{j\omega})| ∣ X ( e j ω ) ∣ 。
根据欧拉公式,正弦信号可展开为两个反向旋转的复指数分量之差:
sin ( ω 0 n ) = e j ω 0 n − e − j ω 0 n 2 j = 1 2 j e j ω 0 n − 1 2 j e − j ω 0 n \sin(\omega_0 n) = \frac{e^{j\omega_0 n} - e^{-j\omega_0 n}}{2j} = \frac{1}{2j} e^{j\omega_0 n} - \frac{1}{2j} e^{-j\omega_0 n} sin ( ω 0 n ) = 2 j e j ω 0 n − e − j ω 0 n = 2 j 1 e j ω 0 n − 2 j 1 e − j ω 0 n
因此,在频域区间 [ − π , π ] [-\pi, \pi] [ − π , π ] 内,幅度谱严格在 ± ω 0 = ± 0.63 rad \pm \omega_0 = \pm 0.63\text{ rad} ± ω 0 = ± 0.63 rad 处隆起两个对称尖峰。
细心的读者会注意到:理论上无限长正弦信号的傅里叶变换应该是两根无限细的冲激 δ \delta δ 函数,为什么图 1 中的峰值呈现出带有起伏波纹的包络?
这正是有限观测带来的第一个物理印记:计算机无法计算负无穷到正无穷的求和,仿真中我们截取了前 32 点数据。截断等价于在理想无限长正弦波上乘以一个矩形窗。时域相乘对应频域卷积,单频脉冲卷积了矩形窗的频谱包络,从而在尖峰两侧拖出了明显的旁瓣。
左下角与右下角分别给出了其相位谱 ∠ X ( e j ω ) \angle X(e^{j\omega}) ∠ X ( e j ω ) 与复平面响应轨迹。由于时域序列的非零支撑区关于原点不对称,截断引入了线性相位倾斜;复平面轨迹则以原点为中心呈现向心回环,直观展示了复数频谱在不同频率下的实部与虚部关联。
2. 自然界的物理约束:实序列的共轭对称性
在工程采集场景中,无论是电极记录的微弱电位、加速度计测得的振动,还是声学传感器捕获的声压,原始信号全部都是实数序列 (x [ n ] ∈ R x[n] \in \mathbb{R} x [ n ] ∈ R )。这个物理前提在数学上必然施加极其严密的对称约束。
我们对 X ( e − j ω ) X(e^{-j\omega}) X ( e − j ω ) 展开推导:
X ( e − j ω ) = ∑ n = − ∞ ∞ x [ n ] e − j ( − ω ) n = ∑ n = − ∞ ∞ x [ n ] e j ω n X(e^{-j\omega}) = \sum_{n=-\infty}^\infty x[n] e^{-j(-\omega)n} = \sum_{n=-\infty}^\infty x[n] e^{j\omega n} X ( e − j ω ) = n = − ∞ ∑ ∞ x [ n ] e − j ( − ω ) n = n = − ∞ ∑ ∞ x [ n ] e j ω n
由于 x [ n ] x[n] x [ n ] 是实数,有 x [ n ] = x ∗ [ n ] x[n] = x^*[n] x [ n ] = x ∗ [ n ] 。对上式取复共轭:
( X ( e − j ω ) ) ∗ = ( ∑ n = − ∞ ∞ x [ n ] e j ω n ) ∗ = ∑ n = − ∞ ∞ x [ n ] e − j ω n = X ( e j ω ) \left( X(e^{-j\omega}) \right)^* = \left( \sum_{n=-\infty}^\infty x[n] e^{j\omega n} \right)^* = \sum_{n=-\infty}^\infty x[n] e^{-j\omega n} = X(e^{j\omega}) ( X ( e − j ω ) ) ∗ = ( n = − ∞ ∑ ∞ x [ n ] e j ω n ) ∗ = n = − ∞ ∑ ∞ x [ n ] e − j ω n = X ( e j ω )
两边再次取共轭,便得到著名的共轭对称性(Conjugate Symmetry) :
X ( e − j ω ) = X ∗ ( e j ω ) X(e^{-j\omega}) = X^*(e^{j\omega}) X ( e − j ω ) = X ∗ ( e j ω )
如果将复数频谱写成极坐标形式 X ( e j ω ) = ∣ X ( e j ω ) ∣ e j θ ( ω ) X(e^{j\omega}) = |X(e^{j\omega})| e^{j \theta(\omega)} X ( e j ω ) = ∣ X ( e j ω ) ∣ e j θ ( ω ) ,代入共轭对称关系可得:
∣ X ( e − j ω ) ∣ = ∣ X ( e j ω ) ∣ ( 幅度偶对称 ) |X(e^{-j\omega})| = |X(e^{j\omega})| \quad (\text{幅度偶对称}) ∣ X ( e − j ω ) ∣ = ∣ X ( e j ω ) ∣ ( 幅度偶对称 )
θ ( − ω ) = − θ ( ω ) ( 相位奇对称 ) \theta(-\omega) = -\theta(\omega) \quad (\text{相位奇对称}) θ ( − ω ) = − θ ( ω ) ( 相位奇对称 )
(图2:实序列频谱的共轭对称性验证)
图 2 给出了双频率叠加实序列的验证结果。在正负对称的数字频率轴 [ − π , π ] [-\pi, \pi] [ − π , π ] 上:
幅度谱以 ω = 0 \omega = 0 ω = 0 为轴严格镜像对称;
相位谱在穿过 ω = 0 \omega = 0 ω = 0 原点时呈现严格的奇对称翻转(正频相位取正,负频相位取负)。
作者按:现实世界中没有虚数信号,负频率并非独立存在的物理实体,而是复指数基底为了抵消虚部、合成纯实数波形而在数学上必须配对的共轭镜像。在高性能计算或嵌入式 DSP 固件开发中,正是由于共轭对称性的存在,实数 FFT(如 CMSIS-DSP 中的 arm_rfft_fast_f32 或科学计算库中的 rfft)只需要计算并存储正频率区间 [ 0 , π ] [0, \pi] [ 0 , π ] 内的 N / 2 + 1 N/2 + 1 N /2 + 1 个复数点,便能完整表达信号信息,直接省下近一半的内存带宽与计算周期。
3. 正交解耦的直观应用:频域滤波与局限性思考
既然信号能在频域被拆解为正交的频率分量,那么信号分离最直接的想法就是:将信号变换到频域,把不需要的频段直接抹掉,再逆变换回时域。
图 3 模拟了这个理想过程。原始信号由 100 Hz 100\text{ Hz} 100 Hz 的基波与 400 Hz 400\text{ Hz} 400 Hz 的高频噪声叠加而成(采样率 f s = 1000 Hz f_s = 1000\text{ Hz} f s = 1000 Hz )。
(图3:理想频域滤波去除高频噪声过程)
在图 3 的右上角,频谱上清晰地竖立着 100 Hz 100\text{ Hz} 100 Hz 与 400 Hz 400\text{ Hz} 400 Hz 两处谱峰。我们构造一个截止频率为 200 Hz 200\text{ Hz} 200 Hz 的理想低通频域掩码 H ( f ) H(f) H ( f ) :在 f ≤ 200 Hz f \le 200\text{ Hz} f ≤ 200 Hz 处保留为 1,在 f > 200 Hz f > 200\text{ Hz} f > 200 Hz 处强制置 0。掩码点乘后,右下角的滤波频谱仅剩下纯净的 100 Hz 100\text{ Hz} 100 Hz 成分;再做逆变换,左下角时域波形恢复为标准的单频正弦波。
但读者诸君在此需要保持严谨的技术警惕:这种直接在频域乘矩形门函数的“理想滤波”,能否直接搬到工程实时数据流中?
答案是否定的。频域的矩形硬截断,其对应的时域冲激响应是一个无限长的 sinc \text{sinc} sinc 函数序列。在实时流式系统中使用有限长块处理强行硬截断,会导致时域波形边界出现严重的吉布斯震荡(Gibbs Phenomenon)和块边缘跳变。这也是为什么我们不能盲目在频域“一刀切”,而必须依赖后续章节中将深入探讨的因果、平滑过渡的 FIR / IIR 滤波器设计。
二、 卷积定理:运算维度的降维打击
在上一篇中,我们推导了线性时不变系统(LTI)的输出等于输入信号 x [ n ] x[n] x [ n ] 与系统单位脉冲响应 h [ n ] h[n] h [ n ] 的卷积和:
y [ n ] = x [ n ] ∗ h [ n ] = ∑ k = − ∞ ∞ x [ k ] h [ n − k ] y[n] = x[n] * h[n] = \sum_{k=-\infty}^\infty x[k] h[n-k] y [ n ] = x [ n ] ∗ h [ n ] = k = − ∞ ∑ ∞ x [ k ] h [ n − k ]
在时域计算该式时,需要进行“翻转、滑动、逐点相乘、求和”四个步骤。如果序列长度较长,计算量随着序列点数急剧增加。
而 DTFT 最核心的性质之一,就是卷积定理(Convolution Property) 。我们给出严密的数学证明:
对输出 y [ n ] y[n] y [ n ] 作 DTFT:
Y ( e j ω ) = ∑ n = − ∞ ∞ y [ n ] e − j ω n = ∑ n = − ∞ ∞ ( ∑ k = − ∞ ∞ x [ k ] h [ n − k ] ) e − j ω n Y(e^{j\omega}) = \sum_{n=-\infty}^{\infty} y[n] e^{-j\omega n} = \sum_{n=-\infty}^{\infty} \left( \sum_{k=-\infty}^{\infty} x[k] h[n-k] \right) e^{-j\omega n} Y ( e j ω ) = n = − ∞ ∑ ∞ y [ n ] e − j ω n = n = − ∞ ∑ ∞ ( k = − ∞ ∑ ∞ x [ k ] h [ n − k ] ) e − j ω n
假设两序列绝对可和,交换双重求和的次序:
Y ( e j ω ) = ∑ k = − ∞ ∞ x [ k ] ( ∑ n = − ∞ ∞ h [ n − k ] e − j ω n ) Y(e^{j\omega}) = \sum_{k=-\infty}^{\infty} x[k] \left( \sum_{n=-\infty}^{\infty} h[n-k] e^{-j\omega n} \right) Y ( e j ω ) = k = − ∞ ∑ ∞ x [ k ] ( n = − ∞ ∑ ∞ h [ n − k ] e − j ω n )
对内层求和引入变元替换,令 m = n − k m = n - k m = n − k (即 n = m + k n = m + k n = m + k ):
∑ n = − ∞ ∞ h [ n − k ] e − j ω n = ∑ m = − ∞ ∞ h [ m ] e − j ω ( m + k ) = e − j ω k ∑ m = − ∞ ∞ h [ m ] e − j ω m = e − j ω k H ( e j ω ) \sum_{n=-\infty}^{\infty} h[n-k] e^{-j\omega n} = \sum_{m=-\infty}^{\infty} h[m] e^{-j\omega (m+k)} = e^{-j\omega k} \sum_{m=-\infty}^{\infty} h[m] e^{-j\omega m} = e^{-j\omega k} H(e^{j\omega}) n = − ∞ ∑ ∞ h [ n − k ] e − j ω n = m = − ∞ ∑ ∞ h [ m ] e − j ω ( m + k ) = e − j ω k m = − ∞ ∑ ∞ h [ m ] e − j ω m = e − j ω k H ( e j ω )
将该项代回外层求和:
Y ( e j ω ) = ∑ k = − ∞ ∞ x [ k ] e − j ω k H ( e j ω ) = ( ∑ k = − ∞ ∞ x [ k ] e − j ω k ) H ( e j ω ) = X ( e j ω ) ⋅ H ( e j ω ) Y(e^{j\omega}) = \sum_{k=-\infty}^{\infty} x[k] e^{-j\omega k} H(e^{j\omega}) = \left( \sum_{k=-\infty}^{\infty} x[k] e^{-j\omega k} \right) H(e^{j\omega}) = X(e^{j\omega}) \cdot H(e^{j\omega}) Y ( e j ω ) = k = − ∞ ∑ ∞ x [ k ] e − j ω k H ( e j ω ) = ( k = − ∞ ∑ ∞ x [ k ] e − j ω k ) H ( e j ω ) = X ( e j ω ) ⋅ H ( e j ω )
由此得到卷积定理的核心等式:
y [ n ] = x [ n ] ∗ h [ n ] ⟺ Y ( e j ω ) = X ( e j ω ) ⋅ H ( e j ω ) y[n] = x[n] * h[n] \quad \Longleftrightarrow \quad Y(e^{j\omega}) = X(e^{j\omega}) \cdot H(e^{j\omega}) y [ n ] = x [ n ] ∗ h [ n ] ⟺ Y ( e j ω ) = X ( e j ω ) ⋅ H ( e j ω )
(图4:卷积定理的时域卷积与频域相乘一致性验证)
图 4 直观呈现了这一数学定理的数值验证:
上排展示了时域计算:输入 x [ n ] = [ 1 , 2 , 3 ] x[n] = [1, 2, 3] x [ n ] = [ 1 , 2 , 3 ] ,系统脉冲响应 h [ n ] = [ 1 , 0.5 ] h[n] = [1, 0.5] h [ n ] = [ 1 , 0.5 ] ,两者的时域卷积输出为 y [ n ] = [ 1 , 2.5 , 4 , 1.5 ] y[n] = [1, 2.5, 4, 1.5] y [ n ] = [ 1 , 2.5 , 4 , 1.5 ] 。
下排展示了频域计算:分别计算 X ( e j ω ) X(e^{j\omega}) X ( e j ω ) 与具有低通平滑特性的 H ( e j ω ) H(e^{j\omega}) H ( e j ω ) ,二者在连续频域逐点相乘得到 ∣ Y ( e j ω ) ∣ |Y(e^{j\omega})| ∣ Y ( e j ω ) ∣ 。
在数值分析中,如果时域序列长度分别为 N N N 和 M M M ,直接时域滑动卷积的计算复杂度为 O ( N ⋅ M ) O(N \cdot M) O ( N ⋅ M ) 。若能借助快速算法转到频域相乘,复杂度将大幅降低。这是频域变换在工业界最重要的计算价值所在。
三、 从连续走向离散的桥梁:引入 DFT 与能量守恒
DTFT 理论极其优美,但它属于数学理想国。它有两个致命缺陷阻碍了其在计算机中的直接落地:
时域积分/求和区间是无限的 (− ∞ -\infty − ∞ 到 + ∞ +\infty + ∞ );
频域变量 ω \omega ω 是连续变量 ,在区间 [ 0 , 2 π ) [0, 2\pi) [ 0 , 2 π ) 内存在不可数无穷多个频率点。
计算机是一台仅能存储有限字长的离散状态机,既无法吃进无限长的数据,也无法画出连续稠密的频谱曲线。
1. 频域采样与时域周期延拓
为了让变换可在计算机中精确计算,工程师采取了两步离散化妥协:
时域截断 :仅截取有限长度为 N N N 的时域序列 x [ n ] x[n] x [ n ] (n = 0 , 1 , … , N − 1 n = 0, 1, \dots, N-1 n = 0 , 1 , … , N − 1 );
频域采样 :在数字角频率的一个完整周期 [ 0 , 2 π ) [0, 2\pi) [ 0 , 2 π ) 内,等间隔均匀抽取 N N N 个离散频点:
ω k = 2 π N k ( k = 0 , 1 , … , N − 1 ) \omega_k = \frac{2\pi}{N} k \quad (k = 0, 1, \dots, N-1) ω k = N 2 π k ( k = 0 , 1 , … , N − 1 )
将 ω k \omega_k ω k 代入 DTFT 定义式,就诞生了离散傅里叶变换(DFT) :
X [ k ] = X ( e j ω ) ∣ ω = 2 π N k = ∑ n = 0 N − 1 x [ n ] e − j 2 π N k n ( k = 0 , 1 , … , N − 1 ) X[k] = X(e^{j\omega}) \Big|_{\omega = \frac{2\pi}{N}k} = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N}kn} \quad (k = 0, 1, \dots, N-1) X [ k ] = X ( e j ω ) ω = N 2 π k = n = 0 ∑ N − 1 x [ n ] e − j N 2 π k n ( k = 0 , 1 , … , N − 1 )
工程中通常引入旋转因子记号 W N = e − j 2 π N W_N = e^{-j\frac{2\pi}{N}} W N = e − j N 2 π ,将其简写为:
X [ k ] = ∑ n = 0 N − 1 x [ n ] W N k n X[k] = \sum_{n=0}^{N-1} x[n] W_N^{kn} X [ k ] = n = 0 ∑ N − 1 x [ n ] W N k n
其对应的 离散傅里叶逆变换(IDFT) 为:
x [ n ] = 1 N ∑ k = 0 N − 1 X [ k ] W N − k n ( n = 0 , 1 , … , N − 1 ) x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] W_N^{-kn} \quad (n = 0, 1, \dots, N-1) x [ n ] = N 1 k = 0 ∑ N − 1 X [ k ] W N − k n ( n = 0 , 1 , … , N − 1 )
(图5:16点矩形脉冲的连续 DTFT 幅度谱与 16 点 DFT 采样点)
图 5 的实验极其深刻:左侧是一个 16 点时域序列(前 4 点为 1,后 12 点为 0 的矩形脉冲);右侧蓝线是其严格由公式计算的连续 DTFT 谱线,而 16 根红色的离散采样针尖正是计算出的 DFT 输出。红色的 DFT 采样点精准地坐落在蓝色的连续 DTFT 曲线上。
但在深入推导之前,必须牢记一个在时域和频域完全对偶的核心法则:
时域连续 ⟺ \Longleftrightarrow ⟺ 频域非周期 ;
时域离散 ⟺ \Longleftrightarrow ⟺ 频域周期延拓 (如 DTFT 频域以 2 π 2\pi 2 π 为周期);
频域离散采样 ⟺ \Longleftrightarrow ⟺ 时域周期延拓 !
当我们在频域以 2 π N \frac{2\pi}{N} N 2 π 为步长等间隔采样出 N N N 个频点 X [ k ] X[k] X [ k ] 时,数学上隐式地对原本有限长的序列 x [ n ] x[n] x [ n ] 施加了周期为 N N N 的周期延拓:
x ~ [ n ] = ∑ r = − ∞ ∞ x [ n − r N ] \tilde{x}[n] = \sum_{r=-\infty}^{\infty} x[n - rN] x ~ [ n ] = r = − ∞ ∑ ∞ x [ n − r N ]
DFT 本质上并不是直接计算孤立的有限长序列,而是计算其周期延拓序列 x ~ [ n ] \tilde{x}[n] x ~ [ n ] 的离散傅里叶级数(DFS)的主值区间!
理解了“频域采样必然诱发时域周期延拓”,你就掌握了破译后文循环卷积、频谱泄漏等一切工程陷阱的总钥匙。
2. 能量不灭定律:帕斯瓦尔定理(Parseval’s Theorem)
时域序列被切断、加权并投影到 N N N 个离散频点后,原信号所蕴含的物理能量是否依然守恒?
设时域序列的总能量定义为各样本模平方和:
E t i m e = ∑ n = 0 N − 1 ∣ x [ n ] ∣ 2 E_{time} = \sum_{n=0}^{N-1} |x[n]|^2 E t im e = n = 0 ∑ N − 1 ∣ x [ n ] ∣ 2
我们将 IDFT 表达式 x [ n ] = 1 N ∑ k = 0 N − 1 X [ k ] W N − k n x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] W_N^{-kn} x [ n ] = N 1 ∑ k = 0 N − 1 X [ k ] W N − k n 代入能量定义式进行严格推导:
∑ n = 0 N − 1 ∣ x [ n ] ∣ 2 = ∑ n = 0 N − 1 x [ n ] x ∗ [ n ] = ∑ n = 0 N − 1 x [ n ] ( 1 N ∑ k = 0 N − 1 X [ k ] W N − k n ) ∗ = 1 N ∑ n = 0 N − 1 x [ n ] ∑ k = 0 N − 1 X ∗ [ k ] W N k n \sum_{n=0}^{N-1} |x[n]|^2 = \sum_{n=0}^{N-1} x[n] x^*[n] = \sum_{n=0}^{N-1} x[n] \left( \frac{1}{N} \sum_{k=0}^{N-1} X[k] W_N^{-kn} \right)^* = \frac{1}{N} \sum_{n=0}^{N-1} x[n] \sum_{k=0}^{N-1} X^*[k] W_N^{kn} n = 0 ∑ N − 1 ∣ x [ n ] ∣ 2 = n = 0 ∑ N − 1 x [ n ] x ∗ [ n ] = n = 0 ∑ N − 1 x [ n ] ( N 1 k = 0 ∑ N − 1 X [ k ] W N − k n ) ∗ = N 1 n = 0 ∑ N − 1 x [ n ] k = 0 ∑ N − 1 X ∗ [ k ] W N k n
交换求和次序:
∑ n = 0 N − 1 ∣ x [ n ] ∣ 2 = 1 N ∑ k = 0 N − 1 X ∗ [ k ] ( ∑ n = 0 N − 1 x [ n ] W N k n ) \sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N} \sum_{k=0}^{N-1} X^*[k] \left( \sum_{n=0}^{N-1} x[n] W_N^{kn} \right) n = 0 ∑ N − 1 ∣ x [ n ] ∣ 2 = N 1 k = 0 ∑ N − 1 X ∗ [ k ] ( n = 0 ∑ N − 1 x [ n ] W N k n )
注意到括号内的求和正是 DFT 的定义 X [ k ] X[k] X [ k ] :
∑ n = 0 N − 1 ∣ x [ n ] ∣ 2 = 1 N ∑ k = 0 N − 1 X ∗ [ k ] X [ k ] = 1 N ∑ k = 0 N − 1 ∣ X [ k ] ∣ 2 \sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N} \sum_{k=0}^{N-1} X^*[k] X[k] = \frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2 n = 0 ∑ N − 1 ∣ x [ n ] ∣ 2 = N 1 k = 0 ∑ N − 1 X ∗ [ k ] X [ k ] = N 1 k = 0 ∑ N − 1 ∣ X [ k ] ∣ 2
这就是离散傅里叶变换的帕斯瓦尔定理 :
∑ n = 0 N − 1 ∣ x [ n ] ∣ 2 = 1 N ∑ k = 0 N − 1 ∣ X [ k ] ∣ 2 \sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2 n = 0 ∑ N − 1 ∣ x [ n ] ∣ 2 = N 1 k = 0 ∑ N − 1 ∣ X [ k ] ∣ 2
(图6:帕斯瓦尔定理数值验证及频域能量分布)
在图 6 的实验中,我们对由三个频率分量叠加的时域信号(左图)进行能量计算:
计算时域总能量:E t i m e = 170.72 E_{time} = 170.72 E t im e = 170.72 ;
计算频域归一化总能量:E f r e q = 1 N ∑ ∣ X [ k ] ∣ 2 = 170.72 E_{freq} = \frac{1}{N}\sum |X[k]|^2 = 170.72 E f r e q = N 1 ∑ ∣ X [ k ] ∣ 2 = 170.72 ;
两者浮点误差为 2.84 × 10 − 14 2.84 \times 10^{-14} 2.84 × 1 0 − 14 ,在 IEEE 754 双精度浮点数的计算舍入边界下严格相等。
为什么公式中频域能量必须除以一个因子 N N N ?
因为 DFT 变换核基底矩阵 W \mathbf{W} W (其中元素为 W N k n W_N^{kn} W N k n )并非酉矩阵(Unitary Matrix)。基底列向量自身的内积为 ∑ n = 0 N − 1 ∣ W N k n ∣ 2 = N \sum_{n=0}^{N-1} |W_N^{kn}|^2 = N ∑ n = 0 N − 1 ∣ W N k n ∣ 2 = N 。这个没有归一化的基底使得变换后的频域矢量模长放大了 N \sqrt{N} N 倍,模平方放大了 N N N 倍。因此只有除以 N N N ,才能真实还原物理能量。
最右侧子图展示的频域能量分布 ∣ X [ k ] ∣ 2 N \frac{|X[k]|^2}{N} N ∣ X [ k ] ∣ 2 ,直接定义了**功率谱密度(PSD)**的离散骨架,是工程中评估各频段信噪比、提取能量特征的基础数学依据。
四、 避坑指南:有限长采样引发的工程幻象
由于 DFT 是对连续频域做等间隔离散采样,并且隐式施加了时域周期延拓,这必然会带来工程上的直觉误区与“视觉幻象”。如果不理解其背后的机理,极易在算法开发与工程调试中写出致命缺陷。
1. 栅栏效应与频谱泄漏:你看到的峰值是真的吗?
当我们用 DFT 观察频谱时,只能看到离散格点 ω k = 2 π N k \omega_k = \frac{2\pi}{N} k ω k = N 2 π k 上的数值。这就像隔着一排栅栏向外张望,视野被固定的一根根木条分割开了,这就是栅栏效应(Picket-fence Effect) 。
如果真实信号的频率恰好落在某个采样格点上,一切看起来风平浪静。但如果真实信号的频率偏离了格点,恰好落在两根栅栏的缝隙中央呢?
图 7 完整揭示了这一场景下的残酷真相。
(图7:格点频率与非格点频率下的栅栏效应及窗函数抑噪表现)
上排:信号频率 f = 78.1 Hz f = 78.1\text{ Hz} f = 78.1 Hz ,恰在格点上 。
采样时长的整数倍恰好包含了信号的完整周期。时域矩形截断在首尾拼接时是平滑的。在右上角的对数频谱中,未加窗(矩形窗)的蓝线呈现出极致干净的单峰,谱线零点恰好穿过所有其他离散格点,看起来“毫无泄漏”。
下排:信号频率 f = 85.9 Hz f = 85.9\text{ Hz} f = 85.9 Hz ,偏离格点半个谱线间距 。
灾难发生了:
峰值被低估(栅栏效应) :由于谱峰最高点落在两根谱线之间,DFT 采到的只是峰尖两侧较低的斜坡,导致读出的主峰幅度被显著削平;
频谱泄漏(Spectral Leakage) :有限长截断导致的 sinc \text{sinc} sinc 旁瓣,其零交叉点不再与 DFT 采样格点对齐。未加窗的矩形窗(蓝线)旁瓣如惊涛骇浪般扩散开来,泄漏电平高达 − 20 dB -20\text{ dB} − 20 dB ,甚至淹没了更远处微弱信号的识别。
窗函数的工程救赎与代价 :
为了压制矩形窗断崖式截断带来的严重旁瓣泄漏,工程上通常在时域引入两端平滑衰减的窗函数(如 Hanning 窗、Hamming 窗、Blackman 窗)。
观察图 7 右下角:
经过加窗后,旁瓣被强力压制到了 − 40 dB -40\text{ dB} − 40 dB 乃至 − 70 dB -70\text{ dB} − 70 dB 以下(红线布莱克曼窗的旁瓣抑制效果最好);
代价是什么?主瓣展宽! 窗函数在时域压低边缘的同时,必然会在频域将其主瓣展宽,降低了分辨两个极度接近频率的能力。
作者按:工程永远是一门权衡的艺术。矩形窗具有最窄的主瓣(最高的频率分辨潜力),却拥有最糟糕的旁瓣(− 13 dB -13\text{ dB} − 13 dB 的高泄漏);Blackman 窗拥有最宽的主瓣(分辨力最差),却能换来极深的高频抑制。在处理包含微弱特征信号的复杂环境时,必须清楚你当前的性能瓶颈究竟是“分不清两个近邻频率”,还是“强信号泄漏掩盖了弱信号”。
2. 补零(Zero-padding)的真相:插值平滑不等于增加分辨率
面对栅栏效应导致的峰值低估,许多初学者经常会采取一种手段:在原始序列末尾补上一大堆 0(Zero-padding) 。
补零后做 FFT,原本稀疏的频域格点瞬间变得密密麻麻,谱线变得极其光滑圆润。有人因此得出结论:“补零提高了频率分辨率”。
在我看来,这是数字信号处理领域流传最广的认知误区之一。
我们必须严格区分两个截然不同的物理概念:
物理频率分辨率(Physical Frequency Resolution) :
Δ f p h y = 1 T t o t a l = f s N r a w \Delta f_{phy} = \frac{1}{T_{total}} = \frac{f_s}{N_{raw}} Δ f p h y = T t o t a l 1 = N r a w f s
它取决于信号的真实有效观测时长 T t o t a l T_{total} T t o t a l 。这是信号本身蕴含的信息量上限。两个频率如果靠得太近,其物理间距小于 1 T t o t a l \frac{1}{T_{total}} T t o t a l 1 ,从信息论的角度看,没有任何后处理算法能够从这一段截断数据中将它们分辨出来。
计算/显示分辨率(Computational / Grid Resolution) :
δ f g r i d = f s N f f t \delta f_{grid} = \frac{f_s}{N_{fft}} δ f g r i d = N f f t f s
通过补零(增大 N f f t N_{fft} N f f t ),我们只是缩小了频域采样的格点间距。
(图8:补零带来的频域插值效应与物理分辨率不变性)
图 8 验证了补零的真实物理作用:
右下图中,原始 32 点信号(蓝线)偏离格点,由于栅栏效应,谱峰读数仅有约 0.35;
当我们在时域尾部补零到 128 点(橙线)和 512 点(绿线)时,高密度的采样点细腻地描绘出了连续 DTFT 的轮廓,精准地找回了被栅栏错过的谱峰幅度(逼近 0.5);
但是请看曲线的宽度 :绿线的主瓣宽度与原始蓝线包络相比,没有丝毫收窄!如果原本有两个距离极近的频率由于有效观测时间过短而融合成一个“大包”,补零后得到的依然只是一个极其平滑的“大包”,绝不可能无中生有地裂变为两个独立的峰。
补零的本质,是在频域做了一次高密度的 sinc 插值,它修复了栅栏效应的采样失真,但绝没有增加任何新的物理信息。
3. 最隐蔽的工程灾难:循环卷积 vs 线性卷积
在实际工业工程中,我们常常需要对两个长序列求时域卷积。根据卷积定理,开发者通常会写出这样的优化逻辑:
时域序列 A, B -> 分别计算 FFT -> 频域逐点相乘 -> 计算 IFFT -> 得到卷积结果
然而,如果不做特殊处理,这段代码输出的结果将引发严重的灾难。因为 DFT 频域相乘逆变换后得到的不是线性卷积,而是循环卷积(Circular Convolution)!
为什么会产生循环卷积?
回顾前文:频域等间隔抽取 N N N 个点,数学上对应将时域序列周期延拓为周期为 N N N 的序列 a ~ [ n ] \tilde{a}[n] a ~ [ n ] 和 b ~ [ n ] \tilde{b}[n] b ~ [ n ] 。
在频域将它们的 N N N 点 DFT 相乘 Y [ k ] = A [ k ] ⋅ B [ k ] Y[k] = A[k] \cdot B[k] Y [ k ] = A [ k ] ⋅ B [ k ] ,在时域对应的正是这两个周期序列在一个周期内的周期卷积 :
y c i r c [ n ] = ∑ m = 0 N − 1 a [ m ] b [ ( n − m ) ( m o d N ) ] ( n = 0 , 1 , … , N − 1 ) y_{circ}[n] = \sum_{m=0}^{N-1} a[m] b[(n - m) \pmod N] \quad (n = 0, 1, \dots, N-1) y c i r c [ n ] = m = 0 ∑ N − 1 a [ m ] b [( n − m ) ( mod N )] ( n = 0 , 1 , … , N − 1 )
如果序列 a [ n ] a[n] a [ n ] 的长度为 N 1 N_1 N 1 ,序列 b [ n ] b[n] b [ n ] 的长度为 N 2 N_2 N 2 ,它们真正的线性卷积 y l i n e a r [ n ] y_{linear}[n] y l in e a r [ n ] 的非零物理支撑区长度为:
L = N 1 + N 2 − 1 L = N_1 + N_2 - 1 L = N 1 + N 2 − 1
当 DFT 点数 N < L N < L N < L 时,由于模 N N N 周期的限制,原本线性卷积中超出前 N N N 个点的部分,会被硬生生回绕折叠(Aliasing)并累加到序列的前部 :
y c i r c [ n ] = ∑ r = − ∞ ∞ y l i n e a r [ n + r N ] ( 0 ≤ n ≤ N − 1 ) y_{circ}[n] = \sum_{r=-\infty}^{\infty} y_{linear}[n + rN] \quad (0 \le n \le N-1) y c i r c [ n ] = r = − ∞ ∑ ∞ y l in e a r [ n + r N ] ( 0 ≤ n ≤ N − 1 )
(图9:循环卷积的时域周期折叠混叠机理与补零法计算线性卷积)
图 9 完整演示了这一混叠碰撞现场:
设 a [ n ] = [ 1 , 2 , 3 , 4 ] a[n] = [1, 2, 3, 4] a [ n ] = [ 1 , 2 , 3 , 4 ] (长 4),b [ n ] = [ 1 , 1 , 1 , 1 ] b[n] = [1, 1, 1, 1] b [ n ] = [ 1 , 1 , 1 , 1 ] (长 4)。
真实的线性卷积长度为 4 + 4 − 1 = 7 4 + 4 - 1 = 7 4 + 4 − 1 = 7 ,理论真实输出为(右上图):
y l i n e a r = [ 1 , 3 , 6 , 10 , 9 , 7 , 4 ] y_{linear} = [1, 3, 6, 10, 9, 7, 4] y l in e a r = [ 1 , 3 , 6 , 10 , 9 , 7 , 4 ]
如果直接用 N = 4 N = 4 N = 4 点计算 DFT 相乘再 IDFT(左下图):
根据回绕折叠法则,超出 N = 4 N=4 N = 4 的第 4、5、6 点(值分别为 9, 7, 4)被折叠累加到第 0、1、2 点上:
y c i r c [ 0 ] = 1 + 9 = 10 y_{circ}[0] = 1 + 9 = 10 y c i r c [ 0 ] = 1 + 9 = 10
y c i r c [ 1 ] = 3 + 7 = 10 y_{circ}[1] = 3 + 7 = 10 y c i r c [ 1 ] = 3 + 7 = 10
y c i r c [ 2 ] = 6 + 4 = 10 y_{circ}[2] = 6 + 4 = 10 y c i r c [ 2 ] = 6 + 4 = 10
y c i r c [ 3 ] = 10 + 0 = 10 y_{circ}[3] = 10 + 0 = 10 y c i r c [ 3 ] = 10 + 0 = 10
计算结果变成了荒谬的全常数序列 [10, 10, 10, 10]!
破解之道 (中下图):
将两个序列在时域末尾补零,使变换长度满足:
N ≥ N 1 + N 2 − 1 N \ge N_1 + N_2 - 1 N ≥ N 1 + N 2 − 1
在本例中,至少补零至 N = 7 N = 7 N = 7 点(工程上为契合 FFT 快速基-2 算法,常补至 ≥ 7 \ge 7 ≥ 7 的 2 的整数次幂,如 8 点)。此时周期的间隔被拉大到了 7 点之外,折叠回绕不再重叠,主值区间内就能丝毫不差地恢复出真正的线性卷积!
五、 物理边界:时宽带宽积与 Dirichlet 核
最后,我们回到时域与频域最深沉的物理羁绊:时间与频率的互斥性 。
设有一个对称有限长矩形脉冲序列:在 − M ≤ n ≤ M -M \le n \le M − M ≤ n ≤ M 时 x [ n ] = 1 x[n] = 1 x [ n ] = 1 ,其余时刻为 0。脉冲的时域总宽度为 2 M + 1 2M + 1 2 M + 1 点。
对其计算 DTFT,推导闭式解:
X ( e j ω ) = ∑ n = − M M e − j ω n = e j ω M ∑ k = 0 2 M e − j ω k = e j ω M 1 − e − j ω ( 2 M + 1 ) 1 − e − j ω X(e^{j\omega}) = \sum_{n=-M}^{M} e^{-j\omega n} = e^{j\omega M} \sum_{k=0}^{2M} e^{-j\omega k} = e^{j\omega M} \frac{1 - e^{-j\omega(2M+1)}}{1 - e^{-j\omega}} X ( e j ω ) = n = − M ∑ M e − j ω n = e j ω M k = 0 ∑ 2 M e − j ω k = e j ω M 1 − e − j ω 1 − e − j ω ( 2 M + 1 )
分子分母分别提取对称相位因子:
X ( e j ω ) = e j ω M e − j ω 2 M + 1 2 ( e j ω 2 M + 1 2 − e − j ω 2 M + 1 2 ) e − j ω 2 ( e j ω 2 − e − j ω 2 ) = sin ( 2 M + 1 2 ω ) sin ( ω 2 ) X(e^{j\omega}) = e^{j\omega M} \frac{e^{-j\omega \frac{2M+1}{2}} \left( e^{j\omega \frac{2M+1}{2}} - e^{-j\omega \frac{2M+1}{2}} \right)}{e^{-j\frac{\omega}{2}} \left( e^{j\frac{\omega}{2}} - e^{-j\frac{\omega}{2}} \right)} = \frac{\sin\left( \frac{2M+1}{2}\omega \right)}{\sin\left( \frac{\omega}{2} \right)} X ( e j ω ) = e j ω M e − j 2 ω ( e j 2 ω − e − j 2 ω ) e − j ω 2 2 M + 1 ( e j ω 2 2 M + 1 − e − j ω 2 2 M + 1 ) = sin ( 2 ω ) sin ( 2 2 M + 1 ω )
这就是信号处理中著名的 Dirichlet 核(周期 sinc 函数) 。
(图10:不同时域宽度矩形脉冲的 DTFT 频谱演化与 Dirichlet 核匹配)
在图 10 中,理论 Dirichlet 曲线(红虚线)与数值计算曲线(蓝实线)严格吻合。我们观察脉冲半宽 M M M 逐渐增大时的演变:
M = 2 M = 2 M = 2 (时域宽 5 点) :主瓣宽达 2.51 rad 2.51\text{ rad} 2.51 rad ,能量在频域广泛扩散;
M = 5 M = 5 M = 5 (时域宽 11 点) :主瓣收窄至 1.14 rad 1.14\text{ rad} 1.14 rad ,主瓣高度上升;
M = 10 M = 10 M = 10 (时域宽 21 点) :主瓣急剧收窄至 0.60 rad 0.60\text{ rad} 0.60 rad ,中心峰值高度攀升至 21。
主瓣的第一零点满足 2 M + 1 2 ω n u l l = π \frac{2M+1}{2}\omega_{null} = \pi 2 2 M + 1 ω n u l l = π ,因此主瓣宽度约为:
Δ ω m a i n = 2 ω n u l l ≈ 4 π 2 M + 1 = 4 π N w i d t h \Delta \omega_{main} = 2\omega_{null} \approx \frac{4\pi}{2M+1} = \frac{4\pi}{N_{width}} Δ ω main = 2 ω n u l l ≈ 2 M + 1 4 π = N w i d t h 4 π
时域有效宽度 N w i d t h N_{width} N w i d t h 与频域主瓣宽度 Δ ω m a i n \Delta \omega_{main} Δ ω main 的乘积是一个固定常数(时宽带宽积常数特性)!
时域上越短暂、越局域化的脉冲,其频域能量越漫溢无边;要想在频域获得一把锐利如手术刀般的窄主瓣探针,唯一的物理途径是在时域维持足够长的物理观测窗口。
这不仅是信号处理的数学定理,更是量子力学海森堡不确定性原理在宏观波动学中的数学同构。
六、 尾声:在物理约束下寻找工程平衡
回顾本篇的探索历程,我们从连续时域的切片开始,借助 DTFT 俯瞰了频域正交投影的本质;通过共轭对称性理解了实数信号的物理约束;用卷积定理建立了时频转化的桥梁;又直面了计算机有限截断与采样所带来的 DFT 妥协,剖析了栅栏效应、频谱泄漏、补零插值和循环卷积混叠这四大工程陷阱。
在我看来,数字信号处理从来不是一门追求“空中楼阁式完美”的纯理论学科,而是一门在严酷物理约束下寻求最优折中的工程艺术。
你无法既要极短的时域响应,又要极窄的频域过渡带;
你无法仅凭后处理的时域补零,就突破物理观测时长的分辨率天花板;
你也无法在不付出计算缓存与延迟的前提下,消除有限字长与截断边界带来的振铃。
承认现实的局限,明晰每一个数学工具背后的假设与代价,并在相互制约的边界中找到最契合业务场景的 Trade-off,这或许就是信号处理工程带给我们的最深层思维训练。
我们下一篇博客见。
附录:本篇可运行验证代码
本篇涉及的全部 10 组仿真实验代码均已开源在代码仓库中:DSP/dsp02/pain_series_02_demo.py 。
# 依赖安装
pip install numpy matplotlib scipy
# 运行本篇仿真实验
python DSP/dsp02/pain_series_02_demo.py