傅里叶变换的深刻思索

2025/12/1910027 字

这篇文章刚好卡住我目前状态的调整阶段,所以真是写了挺久的,不过总算是完成了,希望后面能逐渐恢复写博客的节奏。

有限区间上的无限点函数

最近又看了些傅立叶变换的视频,理解又加深了一些。

可以说,傅立叶变换的本质就是对一个向量进行换基重新表示。

基于这个观点,我们回头看下之前学的原始的傅立叶变换的过程,看能不能从头推导出傅立叶变换的公式。

最开始,我们是对某一段区间之内的函数进行傅立叶变换。

既然我说傅立叶变换就是对向量换基,那么一段区间之内的函数应该如何表示成一个向量呢?

先想想函数的本质是什么,即是什么决定了一个函数是他自己而不是他函数。

函数的本质是从源到像的映射的集合,每一个源只有一个像与之对应,即,只要确定每一对映射,就能确定整个函数。

可以直接把映射的集合作为向量吗?不太行,我们需要建模成具体的数值,才能使用代数的方法进行计算。

那么,由于函数的区间,也就是定义域是确定的,我们只需要把源对应的像的值按照源从大到小的顺序排成一列,也可以确定指定定义域区间上的映射集合,也就确定了函数。

给出如下的数学定义:

函数 f(x)f(x) 在 x∈[a,b]x \in [a, b] 上可以用向量

[f(a)f(a+b−an)⋮f(b)] \begin{bmatrix} f(a) \\ f(a + \frac{b-a}{n}) \\ \vdots \\ f(b) \end{bmatrix}

表示。

当 n→∞n \to \infty 时,向量的长度趋于无穷时,向量的维度无穷,意味着我们取到了函数上的每一点,可以精准地表达一段函数。

好,现在函数可以用一个无限维向量表示了。我们继续用线性代数的方法,把玩一下这个代表函数的向量。

首先是相加相加,发现自然地和函数的加法是匹配的。即两个函数向量相加的结果得到的新向量和两个函数相加得到的函数代表的向量相同:

[f(a)f(a+b−an)⋮f(b)]+[g(a)g(a+b−an)⋮g(b)]=[f(a)+g(a)f(a+b−an)+g(a+b−an)⋮f(b)+g(b)] \begin{bmatrix} f(a) \\ f(a + \frac{b-a}{n}) \\ \vdots \\ f(b) \end{bmatrix} + \begin{bmatrix} g(a) \\ g(a + \frac{b-a}{n}) \\ \vdots \\ g(b) \end{bmatrix} = \begin{bmatrix} f(a) + g(a) \\ f(a + \frac{b-a}{n}) + g(a + \frac{b-a}{n}) \\ \vdots \\ f(b) + g(b) \end{bmatrix} (突然意识到这就是群论中同构的定义,意味着函数向量和函数在加法上是同构的)

继续呢,向量的内积在函数上有何意味?

[f(a)f(a+b−an)⋮f(b)]⋅[g(a)g(a+b−an)⋮g(b)]=∑i=0nf(a+i(b−a)n)g(a+i(b−a)n) \begin{bmatrix} f(a) \\ f(a + \frac{b-a}{n}) \\ \vdots \\ f(b) \end{bmatrix} \cdot \begin{bmatrix} g(a) \\ g(a + \frac{b-a}{n}) \\ \vdots \\ g(b) \end{bmatrix} = \sum_{i=0}^{n} f(a + \frac{i(b-a)}{n}) g(a + \frac{i(b-a)}{n})

当 $ n \to \infty$ 时,右边的式子似乎可以写成积分的形式:

先给右边的式子乘上一个 b−an\frac{b-a}{n}

得到:

lim⁡n→∞b−an∑i=0nf(a+i(b−a)n)g(a+i(b−a)n)=∫abf(x)g(x)dx \lim_{n \to \infty} \frac{b-a}{n} \sum_{i=0}^{n} f(a + \frac{i(b-a)}{n}) g(a + \frac{i(b-a)}{n}) = \int_a^b f(x) g(x) dx

所以:

[f(a)f(a+b−an)⋮f(b)]⋅[g(a)g(a+b−an)⋮g(b)]=nb−a∫abf(x)g(x)dx \begin{bmatrix} f(a) \\ f(a + \frac{b-a}{n}) \\ \vdots \\ f(b) \end{bmatrix} \cdot \begin{bmatrix} g(a) \\ g(a + \frac{b-a}{n}) \\ \vdots \\ g(b) \end{bmatrix} = \frac{n}{b-a} \int_a^b f(x) g(x) dx

结果竟然是两个函数的乘积的积分乘上一个非零常数项,所以,两个函数向量的内积是否为0等价于两个函数乘积在指定区间上的积分是否为0。

这给了我们判断两个函数在向量表示下是否“正交”的方法,因为无限的乘积之和我们算不出来,但是有限区间的积分我们时可以算出来的。

既然我们有办法判断函数对应的无限维向量是否正交了,那为什么不试着找找一组正交基呢,也就是一组向量的集合,其中的向量两两正交。由上述的等价关系,我们的目标就是找到一组函数集合,其中的函数两两之间乘积在指定区间上的积分都是 0。

经过漫长的寻找(其实是傅立叶发现的),我们终于惊喜地发现,还真的存在着这样的一组函数,那就是:

{1,sin(x),cos(x),sin(2x),cos(2x),… } \{1, sin(x), cos(x), sin(2x), cos(2x), \dots\}

具体的证明我就不证了(其实是我不太会手算积分),但是这个结论是正确的,即:

∫02πsin(mx)sin(nx)dx={0m≠nπm=n \int_0^{2\pi} sin(m x) sin(n x) dx = \begin{cases} 0 & m \neq n \\ \pi & m = n \\ \end{cases}
∫02πcos(mx)cos(nx)dx={0m≠nπm=n \int_0^{2\pi} cos(m x) cos(n x) dx = \begin{cases} 0 & m \neq n \\ \pi & m = n\\ \end{cases}
∫02πsin(mx)cos(nx)dx=0 \int_0^{2\pi} sin(m x) cos(n x) dx = 0

那么,有了这样的一组正交基 {1,sin(x),cos(x),sin(2x),cos(2x),… }\{1, sin(x), cos(x), sin(2x), cos(2x), \dots\},我们可不可以试着把某个函数 f(x)f(x) 用这组正交基表示呢?

由格拉姆-施密特正交化公式得:

f(x) 正交化后关于基 g(x)的系数就等于 (用 ⟨a,b⟩\langle a, b \rangle 表示 a,ba, b 做内积运算):

⟨f(x),g(x)⟩⟨g(x),g(x)⟩ \frac{\langle f(x), g(x) \rangle}{\langle g(x), g(x) \rangle}

进一步,使用之前得到的函数向量内积计算的公式:

⟨f(x),g(x)⟩⟨g(x),g(x)⟩=∫02πf(x) g(x) dx∫02πg(x)2dx \frac{\langle f(x), g(x) \rangle}{\langle g(x), g(x) \rangle} = \frac{\int_0^{2\pi} f(x)\, g(x)\, dx}{\int_0^{2\pi} g(x)^2 dx}

这里要注意:不同的基函数 g(x)g(x),分母 ∫02πg(x)2dx\int_0^{2\pi} g(x)^2 dx 不一样。

  • 当 g(x)=1g(x)=1 时,∫02π12dx=2π\int_0^{2\pi} 1^2 dx = 2\pi。
  • 当 g(x)=sin⁡(nx)g(x)=\sin(nx) 或 g(x)=cos⁡(nx)g(x)=\cos(nx)(n≥1n\ge 1)时,∫02πg(x)2dx=π\int_0^{2\pi} g(x)^2 dx = \pi。

因此换基之后的系数(按 {1,sin⁡(x),cos⁡(x),sin⁡(2x),cos⁡(2x),… }\{1,\sin(x),\cos(x),\sin(2x),\cos(2x),\dots\} 排列)应写成:

[12π∫02πf(x)dx1π∫02πf(x)sin(x)dx1π∫02πf(x)cos(x)dx1π∫02πf(x)sin(2x)dx1π∫02πf(x)cos(2x)dx…] \begin{bmatrix} \frac{1}{2\pi}\int_0^{2\pi} f(x) dx \\ \frac{1}{\pi}\int_0^{2\pi} f(x) sin(x) dx \\ \frac{1}{\pi}\int_0^{2\pi} f(x) cos(x) dx \\ \frac{1}{\pi}\int_0^{2\pi} f(x) sin(2x) dx \\ \frac{1}{\pi}\int_0^{2\pi} f(x) cos(2x) dx \\ \dots \end{bmatrix}

有限区间上的离散点函数(或者叫做数列)

沿着之前的思路,还是建模成向量,目前先研究0~2$\pi$区间上的。

这时候突然发现,之前能用的正交基,这个时候不再是正交的了。

比如在[0,2π)[0, 2\pi)区间上,取 4 个点,那就是0,π2,π,3π20, \frac{\pi}{2}, \pi, \frac{3\pi}{2} 这几个点上的函数值。

不难发现 sin(2x)sin(2x) 在这几个点上的取值都是 0,也就是说是 4 维的零向量。这当然不能作为基。

此时再经过一段时间的探索,我们找到了这种情况下的新的基:(应该不是巧合,但是如何从正余弦推到到这里我还不清楚)

{eitx},0≤t≤n−1 \{e^{itx}\}, 0 \le t \le n-1

此时会发现:

⟨eitx,eitx⟩=e2itx \langle e^{itx}, e^{itx} \rangle = e^{2itx}

没用啊,在x=2kπnx = \frac{2k\pi}{n} 上的取值加起来还是为0。

失败了吗?还没有,假如我们把内积定义为(a, b均为n维复数向量):

⟨a,b⟩=∑i=0n−1aibi‾ \langle a, b \rangle = \sum_{i=0}^{n-1} a_i \overline{b_i}

发现此时:

⟨eit1x,eit2x⟩={∑k=0n−1eit12kπneit22kπn‾=∑k=0n−11=nt1=t2∑k=0n−1eit12kπneit22kπn‾=∑k=0n−1eiC2kπn=0t1≠t2 \langle e^{it_1x}, e^{it_2x} \rangle = \begin{cases} \sum_{k=0}^{n-1} e^{it_1\frac{2k\pi}{n}} \overline{e^{it_2\frac{2k\pi}{n}}} = \sum_{k=0}^{n-1} 1 = n & t_1 = t_2 \\ \sum_{k=0}^{n-1} e^{it_1\frac{2k\pi}{n}} \overline{e^{it_2\frac{2k\pi}{n}}} = \sum_{k=0}^{n-1} e^{iC\frac{2k\pi}{n}} = 0 & t_1 \neq t_2 \end{cases}

神奇吧,实数范围上的向量没找到合适的基向量,但是在复数范围上的向量组上找到了。

按照如上的内积的定义,把在[0,2π)在[0, 2\pi)区间上的离散点函数f(x)f(x)换基为{eitx},0≤t≤n−1\{e^{itx}\}, 0 \le t \le n-1,对应的系数依然用格拉姆-施密特方法得到为:

⟨f(x),eitx⟩⟨eitx,eitx⟩=∑k=0n−1f(2kπn)e−it2kπnn,0≤t≤n−1 \frac{\langle f(x), e^{itx} \rangle}{\langle e^{itx}, e^{itx} \rangle} = \frac{\sum_{k=0}^{n-1} f(\frac{2k\pi}{n}) e^{-it\frac{2k\pi}{n}}}{n}, 0 \le t \le n-1

还有更神奇的,让我们把∑k=0n−1f(2kπn)e−it2kπn\sum_{k=0}^{n-1} f(\frac{2k\pi}{n}) e^{-it\frac{2k\pi}{n}}展开:

f(0)+f(2πn)e−it2πn+⋯+f(2(n−1)πn)e−it2(n−1)πn=f(0)(e−it2πn)0+f(2πn)(e−it2πn)1+⋯+f(2(n−1)πn)(e−it2πn)n−1=f(0)x0+f(2πn)x1+⋯+f(2(n−1)πn)xn−1,x=e−it2πn f(0) + f(\frac{2\pi}{n}) e^{-it\frac{2\pi}{n}} + \dots + f(\frac{2(n-1)\pi}{n}) e^{-it\frac{2(n-1)\pi}{n}} \\ = f(0) {(e^{-it\frac{2\pi}{n}})}^0 + f(\frac{2\pi}{n}) {(e^{-it\frac{2\pi}{n}})}^1 + \dots + f(\frac{2(n-1)\pi}{n}) {(e^{-it\frac{2\pi}{n}})}^{n-1} \\ = f(0) x^0 + f(\frac{2\pi}{n}) x^1 + \dots + f(\frac{2(n-1)\pi}{n}) x^{n-1}, x=e^{-it\frac{2\pi}{n}}

看出来了吗?

所以,实际上f(x)f(x)换基后对应的系数就是以数列f(x)f(x)各值为系数的多项式 ∑k=0n−1f(2kπn)xk\sum_{k=0}^{n-1} f(\frac{2k\pi}{n}) x^k 在 x=e−it2πn,0≤t≤n−1x=e^{-it\frac{2\pi}{n}}, 0 \le t \le n-1 的值,再除以n

所以换基的过程可以用一个矩阵表示:

1n[11…11e−i2πn…e−i2(n−1)πn⋮⋮⋱⋮1e−i2(n−1)πn…e−i2(n−1)2πn][f(0)f(2πn)⋮f(2(n−1)πn)]=换基后的系数向量 \frac{1}{n} \begin{bmatrix} 1 & 1 & \dots & 1 \\ 1 & e^{-i\frac{2\pi}{n}} & \dots & e^{-i\frac{2(n-1)\pi}{n}} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & e^{-i\frac{2(n-1)\pi}{n}} & \dots & e^{-i\frac{2(n-1)^2\pi}{n}} \end{bmatrix} \begin{bmatrix} f(0) \\ f(\frac{2\pi}{n}) \\ \vdots \\ f(\frac{2(n-1)\pi}{n}) \end{bmatrix} = 换基后的系数向量

终于可以代码实践一下了,就来画一下4维离散的点的散点图对应的变换的后的函数波形图吧

python
import numpy as np
import matplotlib.pyplot as plt

# 获得 n 维傅立叶变换矩阵
def get_fourier_matrix(n: int):
    return np.array([
        [np.exp(-2j * np.pi * k * t / n) for k in range(n)]
        for t in range(n)
    ])
  
# 获得以e^itx为基,以coeffs为系数的函数
def get_basis_functions(coeffs: np.ndarray):
    def f(x):
        return sum(c * np.exp(1j * t * x) for t, c in enumerate(coeffs))
    return f

# 定义离散点
n = 10
x = [2 * np.pi * k / n for k in range(n)]
y = [np.random.rand() for _ in range(n)]

# 计算傅立叶变换后的系数向量
coeffs = get_fourier_matrix(len(x)) @ y / len(x)
# 获得拟合出来的函数
f = get_basis_functions(coeffs)

# 绘制散点图
plt.scatter(x, y, label='Original')
plt.plot(np.linspace(min(x), max(x), 100), [f(_x) for _x in np.linspace(min(x), max(x), 100)], label='Transformed', color='red')
plt.legend()
plt.show()

图形是对的,说明理论没啥错


快速傅立叶变换(FFT)算法

接下来就是试着优化一下这个变换的算法了,目前我们使用的是矩阵进行变换,也就是说算法的时间复杂度是 O(n2)O(n^2)

还是太慢了,能再快一点吗。

再回顾一下我们把傅立叶变化的任务建模成了什么:

求给定系数的多项式∑k=0n−1f(2kπn)xk\sum_{k=0}^{n-1} f(\frac{2k\pi}{n}) x^k 在 x=e−it2πn,0≤t≤n−1x=e^{-it\frac{2\pi}{n}}, 0 \le t \le n-1 的值,再除以n。

注意到在x=e−it2πn,0≤t≤n−1x=e^{-it\frac{2\pi}{n}}, 0 \le t \le n-1中,有相反数的存在:

即e−i(t+n2)2πn=e−it2πn+π=−e−it2πne^{-i(t+\frac{n}{2})\frac{2\pi}{n}} = e^{-it\frac{2\pi}{n}+\pi} = -e^{-it\frac{2\pi}{n}}

而多项式正好是由偶函数和奇函数构成的,意思是我们只需要算一半的取值就能立即给出全部的取值。

表达成数学公式为,把多项式f(x)分成偶数项和奇数项:

f(x)=feven(x)+fodd(x)f(x) = f_{even}(x) + f_{odd}(x)

又因为 feven(x)=feven(−x)f_{even}(x) = f_{even}(-x),fodd(x)=−fodd(−x)f_{odd}(x) = -f_{odd}(-x),且xi=xi+n2x_i = x_{i+\frac{n}{2}}

得

f(x)=[feven(x[:n2])+fodd(x[:n2])feven(x[:n2])−fodd(x[:n2])] f(x) = \begin{bmatrix} f_{even}(x[:\frac{n}{2}]) + f_{odd}(x[:\frac{n}{2}]) \\ f_{even}(x[:\frac{n}{2}]) - f_{odd}(x[:\frac{n}{2}]) \end{bmatrix}

又注意到 feven(x)=g(x2),fodd(x)=xh(x2)f_{even}(x) = g(x^2), f_{odd}(x) = x h(x^2),其中 g(x2)g(x^2) 是以 f(x)f(x)偶函数项系数为系数的多项式,h(x2)h(x^2) 是以 f(x)f(x)奇函数项系数为系数的多项式。

此时不难注意到我们把原问题拆分成了两个一半规模的子问题。

对子问题g(x2)g(x^2)使用递归,得到的结果是以 f(x)f(x)奇数项系数为系数的多项式在x2=e−it2πn2=e−it4πnx^2=e^{-it\frac{2\pi}{\frac{n}{2}}}=e^{-it\frac{4\pi}{n}},也就是x=e−it2πnx=e^{-it\frac{2\pi}{n}}的值,正好是原问题需要的值。子问题h(x2)h(x^2)也是同理。

所以。递归,启动。

恭喜我们,发现了 FFT(快速傅立叶变换)算法,时间复杂度直接降到 O(nlog(n))O(nlog(n))。

下面是代码实现。

python
import numpy as np
def fft(c: np.ndarray[complex]) -> np.ndarray[complex]:
    if len(c) <= 1:
        return c
    # 补到最近的 2 的幂
    n = 1 << (len(c) - 1).bit_length()
    c = np.concatenate((c, np.zeros(n - len(c), dtype=complex)))

    g = fft(c[::2])
    h = fft(c[1::2])

    res = np.zeros(n, dtype=complex)
    for t in range(n // 2):
        x = np.exp(-2j * np.pi * t / n)
        res[t] = g[t] + x * h[t]
        res[t + n // 2] = g[t] - x * h[t]
    return res

c = np.array([1, 2, 3, 4, 5, 6, 7, 8])
print(np.round(fft(c), 2))

# 比较一下官方的
print(np.round(np.fft.fft(c), 2))

plaintext
[36.+0.j   -4.+9.66j -4.+4.j   -4.+1.66j -4.+0.j   -4.-1.66j -4.-4.j
 -4.-9.66j]
[36.+0.j   -4.+9.66j -4.+4.j   -4.+1.66j -4.+0.j   -4.-1.66j -4.-4.j
 -4.-9.66j]

发现和官方的库求出来答案是一样的,非常好。

逆快速傅立叶变换(IFFT)算法

好,现在反过来,假设我们已经有了傅立叶变换之后得到的系数向量{ct}\{c_t\},我们想要把它变回原来的函数向量。

直接的思路就是,先根据系数向量加上之前定义的基{eitx}\{e^{itx}\},再除以n,得到函数1n∑t=0n−1cieitx\frac{1}{n}\sum_{t=0}^{n-1}{c_ie^{itx}},然后再把函数在指定的离散点{2πkn}\{\frac{2\pi k}{n}\}上取值,得到函数向量。

不过这样的时间复杂度是 O(n2)O(n^2)。变换是快了,但是还原还是原来的速度的话,没啥用。

再思考一下:

∑t=0n−1cieitx\sum_{t=0}^{n-1}{c_ie^{itx}} 实际上也是一个多项式: ∑t=0n−1ciyt,y=eix\sum_{t=0}^{n-1}{c_iy^t}, y=e^{ix}

要求的就是上述的多项式在 x=2πkn,即y=ei2πkn,0≤k≤n−1x=\frac{2\pi k}{n}, 即y=e^{i\frac{2\pi k}{n}}, 0 \le k \le n-1 的值。

同理,也可以使用递归的方式实现 ifft 如下:

python
import numpy as np
def ifft(c: np.ndarray[complex]) -> np.ndarray[complex]:
    if len(c) <= 1:
        return c
    # 补到最近的 2 的幂
    n = 1 << (len(c) - 1).bit_length()
    c = np.concatenate([c, np.zeros(n - len(c), dtype=complex)])

    g = ifft(c[::2])
    h = ifft(c[1::2])

    res = np.zeros(n, dtype=complex)
    for t in range(n // 2):
        y = np.exp(2j * np.pi * t / n) # 只改了这一行
        res[t] = g[t] + y * h[t]
        res[t + n // 2] = g[t] - y * h[t]
    return res

c = np.array([1, 2, 3, 4, 5, 6, 7, 8])
fft_c = fft(c)
ifft_c = ifft(fft_c) / len(fft_c)  # 别忘了除以 n
print(np.round(ifft_c, 2))

# 比较一下官方的
print(np.round(np.fft.ifft(fft_c), 2))
plaintext
[1.+0.j 2.+0.j 3.-0.j 4.+0.j 5.+0.j 6.-0.j 7.+0.j 8.-0.j]
[1.+0.j 2.-0.j 3.-0.j 4.-0.j 5.+0.j 6.+0.j 7.+0.j 8.+0.j]

傅立叶变化的用处

好了,我们明白了傅立叶变换的计算方式,那么,傅立叶变换到底有什么用呢?

再来看一下我们到底把什么东西变成了什么东西

我们把一个在有限区间上的离散点函数,变成了eitxe^{itx}这一组函数的线性组合。

那么,eitxe^{itx}这一组函数有什么特别的地方吗?

注意到,他们都是具有周期T=2πtT=\frac{2\pi}{t}的周期函数,他们的频率不同,但是都是F=tF=t,都是整数。则他们共有的周期是T=2πT=2\pi。

并且,注意到,这一组函数在复平面上的轨迹是一个以原点为中心的圆。

而系数,就是这个圆形轨迹的起点,决定着该圆周运动的振幅和相位。

综上所述,傅立叶变换的结果可以看成是把一个复平面上闭合曲线的轨迹运动,表示成了一组不同频率的圆周运动的叠加。


用傅立叶变换画图

由此而来带来的最直接的一个应用就是用傅立叶变换用首尾相连的圆形运动画出你想要的图形(不过好像也就是娱乐作用了),这部分的作用我这里就不实现了。这部分可以看我以前的博客傅立叶级数


用傅立叶变换进行信号处理

还是回到傅立叶变换的本质上来,我们把一个函数表示成了一组不同频率的圆周运动的叠加。

那么实际上,我们就做到了把一个运动的不同频率的部分进行分离。

在信号处理中,如果我们知道一个信号中哪些频率是我们想要的,哪些频率是噪声,我们就可以通过傅立叶变换把信号分解成不同频率的部分,然后把噪声部分滤掉,再把剩下的部分进行逆傅立叶变换,就能得到一个更干净的信号。

代码实现举例

python
import numpy as np
import matplotlib.pyplot as plt

signal = np.sin(np.linspace(0, 2 * np.pi, 1000))
noise = np.random.normal(loc=0, scale=0.05, size=1000)
noised_signal = signal + noise
fft_signal = np.fft.fft(noised_signal)
filtered_fft_signal = np.array([f if i < 10 else 0 for i, f in enumerate(fft_signal)])
filtered_signal = np.fft.ifft(filtered_fft_signal)

plt.plot(noised_signal)
plt.plot(filtered_signal)
plt.show()

振幅怎么只有一半,怎么回事?

按理来说上述的信号低频部分是绝大多数的啊?而且我们只滤掉后后面990个高频的部分。

那唯一的理由就是,被我们滤掉的信号部分中也有低频的部分。

可是按照上面的数学推导,傅立叶变换后得到的系数就是从0到n-1频率的成分啊,系数的后面怎么会有低频成分呢?

一番推导和查阅资料后,终于发现了原因:

在离散傅立叶变换中,我们用做基的函数是eitx,0≤t≤n−1e^{itx}, 0 \le t \le n-1,此时考虑到x的离散性,x=2πkn,k=0,1,…,n−1x = \frac{2\pi k}{n}, k = 0, 1, \dots, n-1

此时基函数是关于k的数列:eit2πkn,k=0,1,…,n−1e^{it\frac{2\pi k}{n}}, k = 0, 1, \dots, n-1

注意到,当 tt 分别为 t0t_0 和 n−t0n - t_0 时,基函数的取值共轭:

ei(n−t0)2πkn=ei2πk−it02πkn=e−it02πkn e^{i(n - t_0)\frac{2\pi k}{n}} = e^{i2\pi k- i t_0\frac{2\pi k}{n}} = e^{-i t_0 \frac{2\pi k}{n}}

由上述的系数的算法,不难发现,共轭的基函数对应的系数也是共轭的。也就是说,傅立叶变换后计算出的第 t0t_0 个系数和第 n−t0n - t_0 个系数是共轭的。

回到频率的问题。

上述的式子中的e−it2πkn,k=0,1,…,n−1e^{-it\frac{2\pi k}{n}}, k = 0, 1, \dots, n-1从频率的角度也可以理解为,这是一个频率为 −t-t 的复平面上的匀速圆周运动。

直观理解,当 n=4n=4 时,我们有一个频率为 -1 的圆周运动 e−i2πk4,k=0,1,2,3e^{-i\frac{2\pi k}{4}}, k = 0, 1, 2, 3,和一个频率为 3 的圆周运动 ei(3)2πk4,k=0,1,2,3e^{i(3)\frac{2\pi k}{4}}, k = 0, 1, 2, 3,发现这两个离散函数在所有k点上的取值是一样的。也就是说这两个函数描述的运动也是一样的,两者在离散的情况下是等价的。相当于就是频率为 3 的运动一个周期走 3 圈,一个周期有 4 刻,相当于每一刻走 34\frac{3}{4} 圈,而频率为 -1 的运动一个周期走 -1 圈,一个周期有 4 刻,相当于每一刻走 −14-\frac{1}{4} 圈,发现两者每一刻走的圈数其实是一样的。

这下我们就解释了为什么傅立叶变换后得到的系数中,前半部分对应的频率是从 0 到 n2\frac{n}{2},而后半部分对应的频率是从 −n2-\frac{n}{2} 到 -1。也就是说离散傅立叶变换后得到的系数中,前半部分是低频部分,后半部分其实也是低频部分,只不过是负频率而已。

修改代码如下:

python
import numpy as np
import matplotlib.pyplot as plt

n = 1000

signal = np.sin(np.linspace(0, 2 * np.pi, n))
noise = np.random.normal(loc=0, scale=0.1, size=n)
noised_signal = signal + noise
fft_signal = np.fft.fft(noised_signal)
filtered_fft_signal = np.array([f if min(i, n - i) < 10 else 0 for i, f in enumerate(fft_signal)])
filtered_signal = np.fft.ifft(filtered_fft_signal)

plt.plot(noised_signal)
plt.plot(filtered_signal)
plt.show()

这下没问题了

进一步,也可以对二维信号进行傅立叶变换

二维数据最自然的容器就是矩阵:一张灰度图可以写成 f∈Rm×nf \in \mathbb{R}^{m \times n},其中 f(x,y)f(x, y) 表示第 xx 行第 yy 列的像素,再次回到“把函数采样拼成向量”的做法,只不过现在有两条轴需要照顾。

为了沿用“向量空间 + 换基”的套路,先给二维函数定义一个内积。仿照前文的离散情况,直接把矩阵摊平成双重求和:

⟨f,g⟩=∑x=0m−1∑y=0n−1f(x,y) g(x,y)‾ \langle f, g \rangle = \sum_{x=0}^{m-1}\sum_{y=0}^{n-1} f(x, y)\,\overline{g(x, y)}

这个内积满足线性、共轭对称与正定性,因此 (Cm×n,⟨⋅,⋅⟩)(\mathbb{C}^{m \times n}, \langle \cdot, \cdot \rangle) 依旧是希尔伯特空间。接下来照惯例找一组正交基。沿用一维的启发,二维基函数可以写成横纵两个一维复指数的乘积:

φu,v(x,y)=e−i2π(uxm+vyn),0≤u<m,  0≤v<n \varphi_{u,v}(x, y) = e^{-i2\pi (\frac{ux}{m} + \frac{vy}{n})},\quad 0 \le u < m,\; 0 \le v < n

证明这是一组正交基也不难,按照定义:

⟨φu1,v1,φu2,v2⟩=∑x=0m−1∑y=0n−1ei2π(u2−u1)xm+(v2−v1)yn \langle \varphi_{u_1,v_1}, \varphi_{u_2,v_2} \rangle = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1} e^{i2\pi \frac{(u_2-u_1)x}{m} + \frac{(v_2-v_1)y}{n}}

当 (u1,v1)≠(u2,v2)(u_1,v_1) \ne (u_2,v_2) 时,得到的结果为0(易证);当二者完全相等时,两段和都等于各自的维度,结果是 mnmn。因此它们两两正交且张成整个空间,确实构成基。

然后我们来看下这组正交基长什么样子吧

我们知道每一个二维基函数都可以用一个二维矩阵表示,矩阵里面的元素是复数。

那怎么可视化复数呢?这里我们直接取其实部,从 [-1, 1] 映射到 [0, 255] 代表一个灰度像素点

不过要注意的是,上述的基可视化出来的图片只有实部的信息,所以其实是不完整的(也没办法可视化完整),所以不能简单地想象成原图就是上述图片的叠加。

python
import numpy as np
import matplotlib.pyplot as plt
def visualize_basis(m: int, n: int):
    """可视化二维傅立叶基函数"""
    # 先创建一个有 m 行 n 列的画布
    fig, ax = plt.subplots(m, n, figsize=(n, m))
    for u in range(m):
        for v in range(n):
            # 对于每一个基
            result = np.zeros((m, n), dtype=np.uint8)
            for x in range(m):
                for y in range(n):
                    result[x, y] += int((np.real(np.exp(-1j * 2 * np.pi * (u * x / m + v * y / n))) + 1) * 127.5)
  
            # 每个基一张图
            ax[u, v].imshow(result, cmap='gray')
            ax[u, v].axis('off')
    plt.show()

visualize_basis(4, 4)

既然已经有一组正交基,施密特正交化的结果不会改变,只是再次提醒“换到频域坐标的方法是格拉姆-施密特系数公式”:

f^(u,v)=⟨f,φu,v⟩⟨φu,v,φu,v⟩=1mn∑x=0m−1∑y=0n−1f(x,y)e−i2π(uxm+vyn) \hat{f}(u, v) = \frac{\langle f, \varphi_{u,v} \rangle}{\langle \varphi_{u,v}, \varphi_{u,v} \rangle} = \frac{1}{mn} \sum_{x=0}^{m-1}\sum_{y=0}^{n-1} f(x, y) e^{-i2\pi (\frac{ux}{m} + \frac{vy}{n})}

这个表达式就是二维离散傅立叶变换 (2D-DFT)。为了直观感受,用纯 Python 把施密特公式翻译成代码,并和官方 numpy.fft.fft2 对比:

python
import numpy as np

def dft2_manual(data: np.ndarray) -> np.ndarray:
    """直接按施密特公式实现二维离散傅立叶变换"""
    m, n = data.shape
    result = np.zeros((m, n), dtype=np.complex128)
    for u in range(m):  # 枚举频域坐标 u
        for v in range(n):  # 枚举频域坐标 v
            acc = 0.0j
            for x in range(m):  # 在行方向求内积
                for y in range(n):  # 在列方向求内积
                    phase = np.exp(-2j * np.pi * (u * x / m + v * y / n))
                    acc += data[x, y] * phase
            result[u, v] = acc / (m * n)  # 除以基向量的模平方
    return result

demo = np.random.rand(4, 4)
print(np.round(dft2_manual(demo), 2))
print(np.round(np.fft.fft2(demo) / (4 * 4), 2))

plaintext
[[ 0.41+0.j   -0.07-0.05j  0.02-0.j   -0.07+0.05j]
 [ 0.04-0.02j -0.08+0.03j  0.03+0.04j  0.04-0.09j]
 [ 0.02-0.j   -0.03+0.03j -0.05-0.j   -0.03-0.03j]
 [ 0.04+0.02j  0.04+0.09j  0.03-0.04j -0.08-0.03j]]
[[ 0.41+0.j   -0.07-0.05j  0.02+0.j   -0.07+0.05j]
 [ 0.04-0.02j -0.08+0.03j  0.03+0.04j  0.04-0.09j]
 [ 0.02+0.j   -0.03+0.03j -0.05+0.j   -0.03-0.03j]
 [ 0.04+0.02j  0.04+0.09j  0.03-0.04j -0.08-0.03j]]

可以看到自写版本得到的系数与 fft2 结果是一致的,只差一个 mnmn 的尺度因子,这是因为 numpy 默认把 1mn\frac{1}{mn} 留在逆变换中,和我们使用的施密特配方一致。在这一层意义上,两者已经完成了“理念对齐”。

不过就像一维情形那样,双重循环的复杂度是 O(mn⋅mn)O(mn\cdot mn),随便一张图就算到崩溃。优化思路也与一维 FFT 完全平行:

  1. 先把每一行看成独立的一维序列做 FFT。
  2. 再把结果转置/视作按列的序列继续做 FFT。
  3. 这样二维变换的复杂度从 O(N2)O(N^2) 直接降为 O(Nlog⁡N)O(N \log N),其中 N=mnN = mn。

(上面AI写的,着实不容易理解,实际上可以从公式入手)

再看一眼公式,发现

f^(u,v)=1mn∑x=0m−1∑y=0n−1f(x,y)e−i2π(uxm+vyn)=1mn∑x=0m−1e−i2πuxm∑y=0n−1f(x,y)e−i2πvyn \hat{f}(u, v) = \frac{1}{mn} \sum_{x=0}^{m-1}\sum_{y=0}^{n-1} f(x, y) e^{-i2\pi (\frac{ux}{m} + \frac{vy}{n})} = \frac{1}{mn} \sum_{x=0}^{m-1} e^{-i2\pi \frac{ux}{m}} \sum_{y=0}^{n-1} f(x, y) e^{-i2\pi \frac{vy}{n}}

为了和我们代码里的 fft(以及 numpy.fft.fft)对齐,这里把不带归一化的行方向 DFT 记成:

Fx(v)=∑y=0n−1f(x,y)e−i2πvyn F_x(v) = \sum_{y=0}^{n-1} f(x, y) e^{-i2\pi \frac{vy}{n}}

于是上式可以写为:

f^(u,v)=1mn∑x=0m−1e−i2πuxm Fx(v) \hat{f}(u, v) = \frac{1}{mn}\sum_{x=0}^{m-1} e^{-i2\pi \frac{ux}{m}}\, F_x(v)

而 Fx(v)F_x(v) 恰好就是对第 xx 行做一次一维 FFT 后的第 vv 个值:

Fx(v)=fft(f[x,:])[v] F_x(v) = fft(f[x, :])[v]

再令 gv(x)=Fx(v)g_v(x)=F_x(v)(把固定列频率 vv 的结果看成关于 xx 的序列),就得到

∑x=0m−1e−i2πuxm gv(x)=fft([gv(x) for x in range(m)])[u] \sum_{x=0}^{m-1} e^{-i2\pi \frac{ux}{m}}\, g_v(x) = fft([g_v(x)\ for\ x\ in\ range(m)])[u]

因此未归一化的二维变换可以写成:

F(u,v)=∑x=0m−1e−i2πuxm Fx(v)=fft([gv(x) for x in range(m)])[u] F(u,v)=\sum_{x=0}^{m-1} e^{-i2\pi \frac{ux}{m}}\,F_x(v)=fft([g_v(x)\ for\ x\ in\ range(m)])[u]

最后如果要得到前面施密特公式里的 f^(u,v)\hat{f}(u,v),再整体除以 mnmn 即可(即 f^(u,v)=1mnF(u,v)\hat{f}(u,v)=\frac{1}{mn}F(u,v))。

下面的代码直接“借用”我们之前写好的一维 fft/ifft(就在文章前面),然后组合出二维版本

python
import numpy as np
from PIL import Image

def fft2(f: np.ndarray[complex]) -> np.ndarray[complex]:
    """把二维信号矩阵先按行做FFT再按列做FFT"""
    m, n = f.shape
    row_fft = np.array([fft(f[x, :]) for x in range(m)]) # 每一行的fft(这里必须要先算出来存着,不然会很慢)
    col_fft = np.array([fft(row_fft[:, v]) for v in range(n)]).T
    return col_fft

def ifft2(freq: np.ndarray[complex]) -> np.ndarray[complex]:
    """把二维频率矩阵先按行做IFFT再按列做IFFT"""
    m, n = freq.shape
    row_ifft = np.array([ifft(freq[u, :]) for u in range(m)])
    col_ifft = np.array([ifft(row_ifft[:, y]) for y in range(n)]).T
    return col_ifft / (m * n) # 补上除以 mn

# 验证是否和官方库一致
demo = np.array([
    [1, 2, 3, 4],
    [5, 6, 7, 8],
    [9, 10, 11, 12],
    [13, 14, 15, 16]
])
print(f"原矩阵:\n{np.round(demo, 2)}")
print(f"自定义FFT2:\n{np.round(fft2(demo), 2)}")
print(f"官方库FFT2:\n{np.round(np.fft.fft2(demo), 2)}")
print(f"自定义IFFT2:\n{np.round(ifft2(fft2(demo)), 2)}")
print(f"官方库IFFT2:\n{np.round(np.fft.ifft2(fft2(demo)), 2)}")
plaintext
原矩阵:
[[ 1  2  3  4]
 [ 5  6  7  8]
 [ 9 10 11 12]
 [13 14 15 16]]
自定义FFT2:
[[136. +0.j  -8. +8.j  -8. +0.j  -8. -8.j]
 [-32.+32.j   0. +0.j   0. +0.j   0. +0.j]
 [-32. +0.j   0. +0.j   0. +0.j   0. +0.j]
 [-32.-32.j   0. +0.j   0. +0.j   0. +0.j]]
官方库FFT2:
[[136. +0.j  -8. +8.j  -8. +0.j  -8. -8.j]
 [-32.+32.j   0. +0.j   0. +0.j   0. +0.j]
 [-32. +0.j   0. +0.j   0. +0.j   0. +0.j]
 [-32.-32.j   0. +0.j   0. +0.j   0. +0.j]]
自定义IFFT2:
[[ 1.+0.j  2.+0.j  3.+0.j  4.-0.j]
 [ 5.+0.j  6.+0.j  7.+0.j  8.+0.j]
 [ 9.+0.j 10.+0.j 11.+0.j 12.-0.j]
 [13.-0.j 14.-0.j 15.-0.j 16.-0.j]]
官方库IFFT2:
[[ 1.+0.j  2.-0.j  3.+0.j  4.+0.j]
 [ 5.-0.j  6.-0.j  7.-0.j  8.-0.j]
 [ 9.+0.j 10.-0.j 11.+0.j 12.+0.j]
 [13.+0.j 14.+0.j 15.+0.j 16.+0.j]]

最后,用一下自己实现的二维傅立叶变换实现一下图片压缩

python
import numpy as np
from PIL import Image

def compress_image_with_fft2(path: str, keep_ratio: float = 0.1) -> Image.Image:
    """使用自写二维FFT去掉高频块实现图片压缩(不使用 shift 版本)"""
    img = Image.open(path).convert('L')  # 读取灰度图
    data = np.array(img, dtype=np.float32)
    freq = fft2(data)  # 使用自实现二维FFT

    H, W = freq.shape
    # 注意:这里“不做 fftshift”,因此 (0,0) 是直流/低频;矩阵中间区域更接近高频
    # keep_ratio 表示“每个轴上保留的低频比例”(例如 0.1 表示行/列各保留约 10% 的低频带宽)
    ry = max(1, int(H * keep_ratio / 2))
    rx = max(1, int(W * keep_ratio / 2))

    # 保留四个角落的低频块(包含 (0,0) 附近以及环绕对应的“负频”角落),把中间视为高频清零
    filtered = freq.astype(np.complex128, copy=True)
    filtered[ry:H - ry, rx:W - rx] = 0

    recon = ifft2(filtered)  # 逆FFT回到图像域
    recon = np.clip(recon.real, 0, 255).astype(np.uint8)  # 虚部通常是数值误差,取实部并裁剪
    return Image.fromarray(recon)

compressed_image = compress_image_with_fft2('images/lena.png', keep_ratio=0.1).save('images/lena_fft2_compressed.png')

这样从“矩阵视角 → 内积 → 正交基 → 施密特公式 → 分离变量优化 → 图片压缩”整个链条就连起来了:我们先证明二维傅立叶变换依旧是个换基,再用自写 FFT2 把复杂度压下来,最后动手试了下使用fft2进行图片的压缩


用傅立叶变换进行特征提取

终于写到我这次文章的饺子醋了,那就是用傅立叶变换的方式理解卷积神经网络的特征提取,灵感来源于这个视频(参考文献)
卷积神经网络的底层是傅里叶变换,傅里叶变换的底层是希尔伯特空间坐标变换

首先我先说的是,对一个向量进行特征提取,其实也可以看成是对其进行一些变换,使其在变换后的向量空间中相似的向量距离更近,不相似的向量距离更远。从而变得更容易分类。这个变换后的,具有相似度上的语义性的,容易分类的向量,就被叫做特征向量。

那傅立叶变换这种变换能不能给我们提供一种比较好的解释向量的角度呢?

答案是肯定的。

傅立叶变换把一个随时间而变化的在复平面上的坐标函数,表示成了一组不同频率的圆周运动的叠加。让我们可以与时间无关地描述一种运动(因为只需要描述各个不同频率的圆周运动的系数)。也就是说的从时域转到频域。

也就是说,我们可以从原运动(或者说是波形)中提取出与时间无关的特征。

代码举个例子:

python
import numpy as np
import matplotlib.pyplot as plt

def show_data_and_fft_spec(data: np.ndarray):
    fft_data = np.fft.fft(data)
    plt.figure(figsize=(10, 2))
    plt.subplot(1, 2, 1)
    plt.ylim(-0.05, 1.1)
    plt.plot(data)
    plt.title('Original Data')
    plt.subplot(1, 2, 2)
    plt.plot(np.abs(fft_data))
    plt.title('FFT Spectrum')
    plt.show()
    # 保留两位小数
    print(np.round(fft_data, 2))

n = 10

data1 = np.zeros(n)
data1[2:4] = 1
show_data_and_fft_spec(data1)

data2 = np.zeros(n)
data2[6:8] = 1
show_data_and_fft_spec(data2)

data3 = np.zeros(n)
data3[6:8] = 0.5
show_data_and_fft_spec(data3)

data4 = np.zeros(n)
data4[2:4] = 1
data4[6:8] = 0.5
show_data_and_fft_spec(data4)

# 终止plt进程
plt.close('all')

plaintext
[ 2.  +0.j    0.  -1.9j  -1.62+0.j    0.  +1.18j  0.62+0.j    0.  -0.j
  0.62+0.j    0.  -1.18j -1.62+0.j    0.  +1.9j ]

plaintext
[ 2.  +0.j   -1.12+1.54j -0.5 -1.54j  1.12+0.36j -0.5 +0.36j  0.  +0.j
 -0.5 -0.36j  1.12-0.36j -0.5 +1.54j -1.12-1.54j]

plaintext
[ 1.  +0.j   -0.56+0.77j -0.25-0.77j  0.56+0.18j -0.25+0.18j  0.  +0.j
 -0.25-0.18j  0.56-0.18j -0.25+0.77j -0.56-0.77j]

plaintext
[ 3.  +0.j   -0.56-1.13j -1.87-0.77j  0.56+1.36j  0.37+0.18j  0.  -0.j
  0.37-0.18j  0.56-1.36j -1.87+0.77j -0.56+1.13j]

可以看到,虽然 data1 和 data2 在时域上是完全不同的,但是在频域上,他们的傅立叶变换的频谱图是完全一致的。data2 和 data3 在时域上振幅不同,波形相似,在频谱图上也体现了出来。

也就是说,对于一串信号中的一段特定的波形特征,我们有了无视其位置统一的提取方式。

那到底什么地方需要提取波形相似但是和位置无关的特征呢?

答案就是图像的特征提取。

如果我们把图像看成是一个二维的信号,就可以通过傅立叶变换这种类似的方式,提取出图像中的某个特定的纹理特征。

不过,新的问题出现了,在计算傅立叶变换的各个系数时,我们是对整个信号进行计算的,那么如果在信号中不同的地方都出现了一些局部的波形特征(如上述data4),傅立叶变换提取出的波形就不能很好地提取出其中蕴含的多个相似的特征,而是把多个相似特征在不同空间上的拼接判断为一种新的特征。

再举一个上面提到的视频中说的一个例子:这个方式会把数字八的两个圈的组合识别为一种新的特征,而不是两个数字0特征的组合。

那有没有办法解决呢?

有的。


短时傅里叶变换(STFT)

(这块不太懂,让ai沿着我的思路,模仿我的风格写了,感觉还不错)

短时傅里叶变换(STFT)的思路是,把整个信号分成很多个小段,然后对每一小段分别进行傅立叶变换。这样就能提取出每一小段中的局部特征。

还记得前面我们遇到的问题吗?data4 是 data1 和 data3 的组合,但是普通的 FFT 把它识别成了一种全新的特征,而不是两个独立特征的组合。

让我们用 STFT 来解决这个问题:

python
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

def show_stft_comparison(data, title):
    # 计算STFT,窗口重合是为了获得更好的时间分辨率
    f, t, Zxx = signal.stft(data, nperseg=5, noverlap=2)  # 50%重合
  
    plt.figure(figsize=(10, 2))
    plt.subplot(1, 2, 1)
    plt.plot(data)
    plt.title(f'{title} - Time Domain')
    plt.xlabel('Sample')
    plt.ylabel('Amplitude')
  
    plt.subplot(1, 2, 2)
    plt.pcolormesh(t, f, np.abs(Zxx))
    plt.title(f'{title} - STFT Spectrogram')
    plt.xlabel('Time Frame')
    plt.ylabel('Frequency')
    plt.colorbar(label='Magnitude')
  
    plt.tight_layout()
    plt.show()

n = 10

# 重新定义我们之前的测试数据
data1 = np.zeros(n)
data1[2:4] = 1

data3 = np.zeros(n) 
data3[6:8] = 0.5

data4 = np.zeros(n)
data4[2:4] = 1
data4[6:8] = 0.5

# 用STFT分析
show_stft_comparison(data1, 'Data1 (Feature at position 2-3)')
show_stft_comparison(data3, 'Data3 (Feature at position 6-7)')  
show_stft_comparison(data4, 'Data4 (Combined features)')

运行这段代码你会发现一个有趣的现象:

  • data1 的 STFT 热力图在时间轴的前半部分有明显的特征
  • data3 的 STFT 热力图在时间轴的后半部分有明显的特征
  • data4 的 STFT 热力图同时在前半部分和后半部分都有对应的特征!

这就是 STFT 的威力:它能够同时告诉我们"什么特征"在"什么时候"出现 ✨

通过滑动窗口(通常会有重合以获得更好的时间分辨率),STFT 成功地把全局的频域分析变成了局部的时频分析。这样我们就能识别出 data4 实际上是由两个独立的局部特征组成的,而不是一个全新的整体特征。

这个思想后来启发了卷积神经网络中的"卷积核"概念——用小的滤波器在整个输入上滑动,提取局部特征。从数学的角度看,CNN 的卷积操作和 STFT 的滑动窗口在本质上是相通的 🎯

傅立叶变换和卷积的内在联系

首先,卷积是什么?

按照我目前的理解,卷积就是,在两个函数上的运算。具体运算的方式就是在指定的区间内,把一个函数倒过来,和另一个函数进行滑动,把重叠的范围进行相乘再相加的计算。

对于无限个点的连续函数,可以表示为以下的公式:

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

然后卷积关于傅立叶变换有一个卷积定理:

F{f∗g}=F{f}⋅F{g} \mathcal{F}\{f * g\} = \mathcal{F}\{f\} \cdot \mathcal{F}\{g\}

其中 ∗* 表示卷积运算,右边的两个傅立叶变化的结果的点乘表示的是对应的系数直接代数相乘。
卷积定理怎么证我不知道(微积分学的不够好)

注意:在离散情况下,直接使用卷积定理得到的是循环卷积,而不是线性卷积。两个长度为 NN 的数列的线性卷积结果是 2N−12N-1 个点,所以在使用傅立叶变换计算线性卷积时,需要先将两个数列补零到 2N−12N-1 的长度。

具体计算方法:

  1. 补零操作:将原始信号 ff 和 gg 都补零到长度 2N−12N-1

    • 例:f=[1,2,3]f = [1, 2, 3],g=[4,5,6]g = [4, 5, 6]
    • 补零后:f′=[1,2,3,0,0]f' = [1, 2, 3, 0, 0]
    • 补零后:g′=[4,5,6,0,0]g' = [4, 5, 6, 0, 0]
  2. 傅立叶变换:

    • F=FFT(f′)F = \text{FFT}(f'),得到长度为 2N−12N-1 的频谱
    • G=FFT(g′)G = \text{FFT}(g'),得到长度为 2N−12N-1 的频谱
  3. 频域点乘:对 FF 和 GG 的对应系数直接相乘

    • 例:H[k]=F[k]×G[k]H[k] = F[k] \times G[k],k=0,1,...,2N−2k = 0, 1, ..., 2N-2
  4. 逆变换:

    • 对 HH 进行逆傅立叶变换:h=IFFT(H)h = \text{IFFT}(H)
    • 结果 hh 就是线性卷积,长度为 2N−12N-1

感觉傅立叶变换前后的函数在卷积计算上好像是同构的,但是实际论证又不是严格的同构。

然后,从公式上看,卷积的计算和傅立叶变换的计算也非常相似。

总之我想说,傅立叶变换和卷积之间有着非常紧密的联系,而卷积神经网络中的卷积操作可以看作是对输入数据进行局部特征提取的一种方式,这与傅立叶变换在频域上提取信号特征的思想是相通的。

傅立叶变换中我们使用的基函数(如 eitxe^{itx})可以看作是特征提取中的一种卷积核。这也能在一定程度上为我们选择卷积核的大小提供启发——即卷积核的大小要和图像中的基本特征的粒度相匹配,越大的卷积核提取到的特征粒度更大。

💬 评论区

写评论

暂无评论,来发表第一条吧。

© 2026 Corn 的自留地