第 007 篇

傅里叶变换:把一条曲线拆成频率

约 5 分钟 1667 字 信号处理傅里叶变换数学
5 Hz 与 12 Hz 正弦波叠加后的时域波形,以及对应的单边幅度谱
5 Hz 与 12 Hz 正弦波叠加后的时域波形,以及对应的单边幅度谱

一段声音在时间轴上是一条起伏的曲线。同一段声音也可以按频率描述:哪些振动成分比较强,它们的相位如何,彼此叠加后形成什么形状。

傅里叶变换把这两种描述联系起来。它不是只适用于声音的技巧;图像、通信和微分方程里,也能看到同一种“换个基底观察问题”的思路。

为什么会出现复数?

欧拉公式把三角函数装进一个紧凑的表达式:

eiωt=cos⁡(ωt)+isin⁡(ωt).e^{i\omega t}=\cos(\omega t)+i\sin(\omega t).

因此,用复指数表示振动,并没有丢掉正弦与余弦。复数的模与辐角可以分别描述幅度和相位:

z=a+ib,∣z∣=a2+b2,z=∣z∣eiarg⁡z.z=a+ib,\qquad |z|=\sqrt{a^2+b^2},\qquad z=|z|e^{i\arg z}.

如果只保留频谱的模,通常就无法唯一还原原信号。两段信号可以有相同的幅度谱,却因为相位不同而呈现不同的波形。

先固定一套变换约定

本文采用角频率 ω\omega,它与通常以赫兹表示的频率 ν\nu 满足:

ω=2πν.\omega=2\pi\nu.

定义傅里叶变换为:

f^(ω)=∫−∞∞f(t)e−iωt dt.\boxed{\widehat f(\omega)= \int_{-\infty}^{\infty}f(t)e^{-i\omega t}\,dt}.

在满足相应反演条件时,逆变换为:

f(t)=12π∫−∞∞f^(ω)eiωt dω.\boxed{f(t)=\frac1{2\pi} \int_{-\infty}^{\infty}\widehat f(\omega)e^{i\omega t}\,d\omega}.

例如,ff 与 f^\widehat f 都可积时,可以在适当意义下使用反演公式;对平方可积信号,则可以通过 L2L^2 理论扩展。后面提到的导数等性质,也各自需要相应的正则性和边界条件。

不同书籍可能把 2π2\pi 分配在不同位置。比较公式前,先检查频率变量与归一化约定。[1]

一个矩形脉冲,会变成什么?

考虑宽度为 T>0T>0、高度为一的矩形脉冲:

f(t)={1,∣t∣≤T/2,0,∣t∣>T/2.f(t)= \begin{cases} 1,& |t|\le T/2,\\ 0,& |t|>T/2. \end{cases}

代入定义:

f^(ω)=∫−T/2T/2e−iωt dt=2sin⁡(ωT/2)ω=T sinc⁡ ⁣(ωT2).\begin{aligned} \widehat f(\omega) &=\int_{-T/2}^{T/2}e^{-i\omega t}\,dt\\ &=\frac{2\sin(\omega T/2)}{\omega}\\ &=T\,\operatorname{sinc}\!\left(\frac{\omega T}{2}\right). \end{aligned}

这里使用的定义是 sinc⁡(x)=sin⁡x/x\operatorname{sinc}(x)=\sin x/x,并通过连续延拓规定 sinc⁡(0)=1\operatorname{sinc}(0)=1,所以 f^(0)=T\widehat f(0)=T。

脉冲在时间上截得很整齐,频谱却会向外延伸。越窄的时间脉冲,会对应越宽的频谱。这也解释了为什么“只看有限一段数据”会影响频率分析。

一个更平滑的例子是高斯函数:

f(t)=e−t2/(2σ2)⟹f^(ω)=σ2π e−σ2ω2/2,σ>0.f(t)=e^{-t^2/(2\sigma^2)} \quad\Longrightarrow\quad \widehat f(\omega)= \sigma\sqrt{2\pi}\,e^{-\sigma^2\omega^2/2}, \qquad \sigma>0.

它变换以后仍是高斯形状,时间宽度与频率宽度则朝相反方向变化。

几条性质,能省掉大量计算

傅里叶变换具有线性:

F{af+bg}=af^+bg^.\mathcal F\{af+bg\} =a\widehat f+b\widehat g.

把信号延迟 t0t_0,频域里会出现一个相位因子:

F{f(t−t0)}=e−iωt0f^(ω).\mathcal F\{f(t-t_0)\} =e^{-i\omega t_0}\widehat f(\omega).

在函数足够光滑且边界项消失时,求导对应乘以 iωi\omega:

F{f′(t)}=iωf^(ω).\mathcal F\{f'(t)\}=i\omega\widehat f(\omega).

定义卷积:

(f∗g)(t)=∫−∞∞f(τ)g(t−τ) dτ,(f*g)(t)=\int_{-\infty}^{\infty}f(\tau)g(t-\tau)\,d\tau,

那么,在卷积定理适用的条件下:

F{f∗g}=f^ g^.\mathcal F\{f*g\}=\widehat f\,\widehat g.

这使得滤波很容易理解:时间域里与某个核做卷积,相当于频率域里对不同成分乘以不同权重。

对应本文的归一化约定,能量关系是:

∫−∞∞∣f(t)∣2 dt=12π∫−∞∞∣f^(ω)∣2 dω.\int_{-\infty}^{\infty}|f(t)|^2\,dt = \frac1{2\pi} \int_{-\infty}^{\infty}|\widehat f(\omega)|^2\,d\omega.

它说明转换描述方式并没有凭空创造或消灭能量。

电脑拿到的是样本

连续积分适合理论分析。电脑实际得到的通常是 NN 个样本:

xn=f(nΔt),n=0,1,…,N−1,fs=1Δt.x_n=f(n\Delta t),\qquad n=0,1,\ldots,N-1,\qquad f_s=\frac1{\Delta t}.

离散傅里叶变换,也就是 DFT,定义为:

Xk=∑n=0N−1xne−i2πkn/N,k=0,…,N−1.X_k=\sum_{n=0}^{N-1}x_n e^{-i2\pi kn/N}, \qquad k=0,\ldots,N-1.

它的逆变换为:

xn=1N∑k=0N−1Xkei2πkn/N.x_n=\frac1N\sum_{k=0}^{N-1}X_k e^{i2\pi kn/N}.

正向与逆向的符号和 1/N1/N 要配套。常见的 FFT 是快速计算 DFT 的算法,并不是另一种变换。[2]

把 DFT 看成有限时间区间内连续积分的近似时,还会涉及 Δt\Delta t 的比例因子;直接把 DFT 数组当成前面连续变换的数值,容易弄错单位和幅度。

两个频率,做一个可复现的例子

构造一个合成信号:

x(t)=sin⁡(2π⋅5t)+0.5sin⁡(2π⋅12t).x(t)=\sin(2\pi\cdot5t)+0.5\sin(2\pi\cdot12t).

它由幅度为一的 5 Hz5\,\mathrm{Hz} 正弦波,与幅度为 0.50.5 的 12 Hz12\,\mathrm{Hz} 正弦波相加。使用 fs=128 Hzf_s=128\,\mathrm{Hz},采集 N=128N=128 个样本:

Δν=fsN=1 Hz.\Delta\nu=\frac{f_s}{N}=1\,\mathrm{Hz}.

这两个频率都恰好落在 DFT 的频点上,适合观察幅度标定。对于实信号、偶数 NN,单边幅度谱的内部频点可计算为:

Ak=2∣Xk∣N,1≤k<N/2.A_k=\frac{2|X_k|}{N}, \qquad 1\le k<N/2.

直流分量 k=0k=0 和奈奎斯特频点 k=N/2k=N/2 不乘二。

import numpy as np

fs = 128
n = 128
t = np.arange(n) / fs
x = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 12 * t)

spectrum = np.fft.rfft(x)
frequency = np.fft.rfftfreq(n, d=1 / fs)
amplitude = np.abs(spectrum) / n
amplitude[1:-1] *= 2  # 此例 n 为偶数;不加倍直流与奈奎斯特频点

print(frequency[5], round(amplitude[5], 6))    # 5.0 1.0
print(frequency[12], round(amplitude[12], 6)) # 12.0 0.5

上图是两个正弦分量叠加后的时域波形,下图在 5 Hz 和 12 Hz 处分别出现幅度 1 和 0.5 的谱峰

这张图由同样的合成信号生成;对应的两个谱峰已通过数值计算核对。它不是现场录音的频谱。

最常见的两个误会

采样够快,还需要知道信号范围

对理想带限信号,通常要求采样频率严格高于最高频率的两倍:

fs>2νmax⁡.f_s>2\nu_{\max}.

若高频成分在采样前没有被限制,它们可能混叠成低频。采样之后再提高 FFT 点数,无法恢复已经混叠的信息。

频谱旁边的波纹,不一定是新频率

有限长度的数据,相当于给原信号乘了一个时间窗。在本文约定下,时间域相乘对应频域卷积:

F{fw}=12π(f^∗w^).\mathcal F\{fw\} =\frac1{2\pi}(\widehat f*\widehat w).

当观测区间不能容纳整数个周期时,能量可能散到相邻频点,形成频谱泄漏。加窗可以降低某些旁瓣,但往往会加宽主瓣,也会影响幅度,需要按窗函数做相应校正。

补零能够让频谱曲线采样得更密,却不会增加观测时间,也不会自动提高分辨两个邻近频率的能力。

理解傅里叶变换之后,再看一条起伏曲线,就多了一种问题可问:它由哪些尺度的变化组成?这些变化是信号自身的结构,还是观测方式留下的痕迹?频谱有用,前提是我们也记得它怎样被算出来。

延伸阅读

  1. Fourier transform
  2. NumPy:离散傅里叶变换与归一化约定

评论

还没有评论。 来写第一条。

留下一条评论

选填,不会公开

纯文本,链接会自动识别。