傅里叶变换的深刻思索
这篇文章刚好卡住我目前状态的调整阶段,所以真是写了挺久的,不过总算是完成了,希望后面能逐渐恢复写博客的节奏。
有限区间上的无限点函数
最近又看了些傅立叶变换的视频,理解又加深了一些。
可以说,傅立叶变换的本质就是对一个向量进行换基重新表示。
基于这个观点,我们回头看下之前学的原始的傅立叶变换的过程,看能不能从头推导出傅立叶变换的公式。
最开始,我们是对某一段区间之内的函数进行傅立叶变换。
既然我说傅立叶变换就是对向量换基,那么一段区间之内的函数应该如何表示成一个向量呢?
先想想函数的本质是什么,即是什么决定了一个函数是他自己而不是他函数。
函数的本质是从源到像的映射的集合,每一个源只有一个像与之对应,即,只要确定每一对映射,就能确定整个函数。
可以直接把映射的集合作为向量吗?不太行,我们需要建模成具体的数值,才能使用代数的方法进行计算。
那么,由于函数的区间,也就是定义域是确定的,我们只需要把源对应的像的值按照源从大到小的顺序排成一列,也可以确定指定定义域区间上的映射集合,也就确定了函数。
给出如下的数学定义:
函数
表示。
当
好,现在函数可以用一个无限维向量表示了。我们继续用线性代数的方法,把玩一下这个代表函数的向量。
首先是相加相加,发现自然地和函数的加法是匹配的。即两个函数向量相加的结果得到的新向量和两个函数相加得到的函数代表的向量相同:
继续呢,向量的内积在函数上有何意味?
当 $ n \to \infty$ 时,右边的式子似乎可以写成积分的形式:
先给右边的式子乘上一个
得到:
所以:
结果竟然是两个函数的乘积的积分乘上一个非零常数项,所以,两个函数向量的内积是否为0等价于两个函数乘积在指定区间上的积分是否为0。
这给了我们判断两个函数在向量表示下是否“正交”的方法,因为无限的乘积之和我们算不出来,但是有限区间的积分我们时可以算出来的。
既然我们有办法判断函数对应的无限维向量是否正交了,那为什么不试着找找一组正交基呢,也就是一组向量的集合,其中的向量两两正交。由上述的等价关系,我们的目标就是找到一组函数集合,其中的函数两两之间乘积在指定区间上的积分都是 0。
经过漫长的寻找(其实是傅立叶发现的),我们终于惊喜地发现,还真的存在着这样的一组函数,那就是:
具体的证明我就不证了(其实是我不太会手算积分),但是这个结论是正确的,即:
那么,有了这样的一组正交基
由格拉姆-施密特正交化公式得:
f(x) 正交化后关于基 g(x)的系数就等于 (用
进一步,使用之前得到的函数向量内积计算的公式:
这里要注意:不同的基函数
- 当
时, 。 - 当
或 ( )时, 。
因此换基之后的系数(按
有限区间上的离散点函数(或者叫做数列)
沿着之前的思路,还是建模成向量,目前先研究0~2$\pi$区间上的。
这时候突然发现,之前能用的正交基,这个时候不再是正交的了。
比如在
不难发现
此时再经过一段时间的探索,我们找到了这种情况下的新的基:(应该不是巧合,但是如何从正余弦推到到这里我还不清楚)
此时会发现:
没用啊,在
失败了吗?还没有,假如我们把内积定义为(a, b均为n维复数向量):
发现此时:
神奇吧,实数范围上的向量没找到合适的基向量,但是在复数范围上的向量组上找到了。
按照如上的内积的定义,把
还有更神奇的,让我们把
看出来了吗?
所以,实际上
所以换基的过程可以用一个矩阵表示:
终于可以代码实践一下了,就来画一下4维离散的点的散点图对应的变换的后的函数波形图吧
pythonimport 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)算法
接下来就是试着优化一下这个变换的算法了,目前我们使用的是矩阵进行变换,也就是说算法的时间复杂度是
还是太慢了,能再快一点吗。
再回顾一下我们把傅立叶变化的任务建模成了什么:
求给定系数的多项式
注意到在
即
而多项式正好是由偶函数和奇函数构成的,意思是我们只需要算一半的取值就能立即给出全部的取值。
表达成数学公式为,把多项式f(x)分成偶数项和奇数项:
又因为
得
又注意到
此时不难注意到我们把原问题拆分成了两个一半规模的子问题。
对子问题
所以。递归,启动。
恭喜我们,发现了 FFT(快速傅立叶变换)算法,时间复杂度直接降到
下面是代码实现。
pythonimport 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)算法
好,现在反过来,假设我们已经有了傅立叶变换之后得到的系数向量
直接的思路就是,先根据系数向量加上之前定义的基
不过这样的时间复杂度是
再思考一下:
要求的就是上述的多项式在
同理,也可以使用递归的方式实现 ifft 如下:
pythonimport 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]
傅立叶变化的用处
好了,我们明白了傅立叶变换的计算方式,那么,傅立叶变换到底有什么用呢?
再来看一下我们到底把什么东西变成了什么东西
我们把一个在有限区间上的离散点函数,变成了
那么,
注意到,他们都是具有周期
并且,注意到,这一组函数在复平面上的轨迹是一个以原点为中心的圆。
而系数,就是这个圆形轨迹的起点,决定着该圆周运动的振幅和相位。
综上所述,傅立叶变换的结果可以看成是把一个复平面上闭合曲线的轨迹运动,表示成了一组不同频率的圆周运动的叠加。
用傅立叶变换画图
由此而来带来的最直接的一个应用就是用傅立叶变换用首尾相连的圆形运动画出你想要的图形(不过好像也就是娱乐作用了),这部分的作用我这里就不实现了。这部分可以看我以前的博客傅立叶级数
用傅立叶变换进行信号处理
还是回到傅立叶变换的本质上来,我们把一个函数表示成了一组不同频率的圆周运动的叠加。
那么实际上,我们就做到了把一个运动的不同频率的部分进行分离。
在信号处理中,如果我们知道一个信号中哪些频率是我们想要的,哪些频率是噪声,我们就可以通过傅立叶变换把信号分解成不同频率的部分,然后把噪声部分滤掉,再把剩下的部分进行逆傅立叶变换,就能得到一个更干净的信号。
代码实现举例
pythonimport 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频率的成分啊,系数的后面怎么会有低频成分呢?
一番推导和查阅资料后,终于发现了原因:
在离散傅立叶变换中,我们用做基的函数是
此时基函数是关于k的数列:
注意到,当
由上述的系数的算法,不难发现,共轭的基函数对应的系数也是共轭的。也就是说,傅立叶变换后计算出的第
回到频率的问题。
上述的式子中的
直观理解,当
这下我们就解释了为什么傅立叶变换后得到的系数中,前半部分对应的频率是从 0 到
修改代码如下:
pythonimport 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()
这下没问题了
进一步,也可以对二维信号进行傅立叶变换
二维数据最自然的容器就是矩阵:一张灰度图可以写成
为了沿用“向量空间 + 换基”的套路,先给二维函数定义一个内积。仿照前文的离散情况,直接把矩阵摊平成双重求和:
这个内积满足线性、共轭对称与正定性,因此
证明这是一组正交基也不难,按照定义:
当
然后我们来看下这组正交基长什么样子吧
我们知道每一个二维基函数都可以用一个二维矩阵表示,矩阵里面的元素是复数。
那怎么可视化复数呢?这里我们直接取其实部,从 [-1, 1] 映射到 [0, 255] 代表一个灰度像素点
不过要注意的是,上述的基可视化出来的图片只有实部的信息,所以其实是不完整的(也没办法可视化完整),所以不能简单地想象成原图就是上述图片的叠加。
pythonimport 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)
既然已经有一组正交基,施密特正交化的结果不会改变,只是再次提醒“换到频域坐标的方法是格拉姆-施密特系数公式”:
这个表达式就是二维离散傅立叶变换 (2D-DFT)。为了直观感受,用纯 Python 把施密特公式翻译成代码,并和官方 numpy.fft.fft2 对比:
pythonimport 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 结果是一致的,只差一个
不过就像一维情形那样,双重循环的复杂度是
- 先把每一行看成独立的一维序列做 FFT。
- 再把结果转置/视作按列的序列继续做 FFT。
- 这样二维变换的复杂度从
直接降为 ,其中 。
(上面AI写的,着实不容易理解,实际上可以从公式入手)
再看一眼公式,发现
为了和我们代码里的 fft(以及 numpy.fft.fft)对齐,这里把不带归一化的行方向 DFT 记成:
于是上式可以写为:
而
再令
因此未归一化的二维变换可以写成:
最后如果要得到前面施密特公式里的
下面的代码直接“借用”我们之前写好的一维 fft/ifft(就在文章前面),然后组合出二维版本
pythonimport 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]]
最后,用一下自己实现的二维傅立叶变换实现一下图片压缩
pythonimport 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进行图片的压缩
用傅立叶变换进行特征提取
终于写到我这次文章的饺子醋了,那就是用傅立叶变换的方式理解卷积神经网络的特征提取,灵感来源于这个视频(参考文献)
卷积神经网络的底层是傅里叶变换,傅里叶变换的底层是希尔伯特空间坐标变换
首先我先说的是,对一个向量进行特征提取,其实也可以看成是对其进行一些变换,使其在变换后的向量空间中相似的向量距离更近,不相似的向量距离更远。从而变得更容易分类。这个变换后的,具有相似度上的语义性的,容易分类的向量,就被叫做特征向量。
那傅立叶变换这种变换能不能给我们提供一种比较好的解释向量的角度呢?
答案是肯定的。
傅立叶变换把一个随时间而变化的在复平面上的坐标函数,表示成了一组不同频率的圆周运动的叠加。让我们可以与时间无关地描述一种运动(因为只需要描述各个不同频率的圆周运动的系数)。也就是说的从时域转到频域。
也就是说,我们可以从原运动(或者说是波形)中提取出与时间无关的特征。
代码举个例子:
pythonimport 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 来解决这个问题:
pythonimport 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 的滑动窗口在本质上是相通的 🎯
傅立叶变换和卷积的内在联系
首先,卷积是什么?
按照我目前的理解,卷积就是,在两个函数上的运算。具体运算的方式就是在指定的区间内,把一个函数倒过来,和另一个函数进行滑动,把重叠的范围进行相乘再相加的计算。
对于无限个点的连续函数,可以表示为以下的公式:
然后卷积关于傅立叶变换有一个卷积定理:
其中
卷积定理怎么证我不知道(微积分学的不够好)
注意:在离散情况下,直接使用卷积定理得到的是循环卷积,而不是线性卷积。两个长度为
具体计算方法:
-
补零操作:将原始信号
和 都补零到长度 - 例:
, - 补零后:
- 补零后:
- 例:
-
傅立叶变换:
,得到长度为 的频谱 ,得到长度为 的频谱
-
频域点乘:对
和 的对应系数直接相乘 - 例:
,
- 例:
-
逆变换:
- 对
进行逆傅立叶变换: - 结果
就是线性卷积,长度为
- 对
感觉傅立叶变换前后的函数在卷积计算上好像是同构的,但是实际论证又不是严格的同构。
然后,从公式上看,卷积的计算和傅立叶变换的计算也非常相似。
总之我想说,傅立叶变换和卷积之间有着非常紧密的联系,而卷积神经网络中的卷积操作可以看作是对输入数据进行局部特征提取的一种方式,这与傅立叶变换在频域上提取信号特征的思想是相通的。
傅立叶变换中我们使用的基函数(如
💬 评论区
暂无评论,来发表第一条吧。