别怕数学！用Python从零实现图像傅里叶变换（附完整代码与频谱图分析）

当你第一次听说“傅里叶变换”时，脑海中是不是立刻浮现出一堆令人望而生畏的数学符号？别担心，这篇文章将用最直观的方式带你理解这个看似高深的概念。就像厨师通过品尝就能分辨出一道菜的配料比例，傅里叶变换让我们能够“品尝”图像的频率成分——这正是图像处理的魔法所在。

我们将完全从零开始，用Python和NumPy一步步构建自己的傅里叶变换实现。没有复杂的公式推导，只有清晰的代码和可视化的频谱图。读完本文后，你将能够：

• 理解傅里叶变换如何将图像分解为不同频率的波
• 手动实现一维和二维离散傅里叶变换（DFT）
• 使用频谱图直观分析图像特征
• 比较自实现与NumPy内置FFT的性能差异

1. 傅里叶变换的直观理解

想象你正在听一首交响乐。虽然听到的是整体声音，但实际上它由许多不同频率的乐器声组合而成。傅里叶变换就像是一个“音乐分解器”，能够将复杂的交响乐分解为各个乐器的单独频率成分。

对于图像处理也是如此。任何图像都可以看作是由不同频率的“波”组成的：

• 低频成分：对应图像中缓慢变化的区域，如大面积的天空或墙面
• 高频成分：对应图像中快速变化的边缘和细节

1.1 从正弦波到图像分解

让我们用简单的正弦波来建立直觉。下面这段代码生成三个不同频率的正弦波，并将它们相加：

import numpy as np
import matplotlib.pyplot as plt

# 生成三个不同频率的正弦波
x = np.linspace(0, 2*np.pi, 1000)
wave1 = 0.5 * np.sin(1 * x)  # 低频
wave2 = 0.3 * np.sin(5 * x)  # 中频
wave3 = 0.2 * np.sin(10 * x) # 高频

# 组合波
combined = wave1 + wave2 + wave3

# 可视化
plt.figure(figsize=(12, 8))
plt.subplot(4, 1, 1); plt.plot(x, wave1); plt.title('低频波')
plt.subplot(4, 1, 2); plt.plot(x, wave2); plt.title('中频波')
plt.subplot(4, 1, 3); plt.plot(x, wave3); plt.title('高频波')
plt.subplot(4, 1, 4); plt.plot(x, combined); plt.title('组合波')
plt.tight_layout()
plt.show()

傅里叶变换的逆过程就是：给定组合波，如何找出其中的wave1、wave2和wave3？这正是我们要实现的核心功能。

2. 一维离散傅里叶变换实现

现在让我们动手实现一维DFT。虽然NumPy已经有现成的fft函数，但自己实现一次会让你真正理解其工作原理。

2.1 DFT的矩阵形式

离散傅里叶变换可以用矩阵乘法表示：

F = W · f

其中：

• f是输入信号（长度为N的向量）
• W是N×N的变换矩阵
• F是输出的频域表示

变换矩阵W的每个元素为：

W[m, n] = e^(-2πi * m*n / N)

这里i是虚数单位。用Python实现如下：

def dft_1d(signal):
    N = len(signal)
    n = np.arange(N)
    k = n.reshape((N, 1))
    
    # 创建变换矩阵
    W = np.exp(-2j * np.pi * k * n / N)
    
    # 矩阵乘法得到频域表示
    freq_domain = np.dot(W, signal)
    return freq_domain

2.2 测试一维DFT

让我们用这个函数分析一个简单信号：

# 创建测试信号：两个正弦波的组合
N = 1000
t = np.linspace(0, 1, N)
signal = 0.5*np.sin(2*np.pi*5*t) + 0.3*np.sin(2*np.pi*20*t)

# 计算DFT
freq_domain = dft_1d(signal)

# 计算频率轴
freqs = np.fft.fftfreq(N, d=1/N)

# 可视化
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(t, signal)
plt.title('时域信号')

plt.subplot(1, 2, 2)
plt.plot(freqs[:N//2], np.abs(freq_domain[:N//2]))
plt.title('频域表示')
plt.xlabel('频率(Hz)')
plt.tight_layout()
plt.show()

你会看到频域图中在5Hz和20Hz处有两个明显的峰值，这正是我们合成信号中使用的频率。

3. 二维图像傅里叶变换

图像是二维信号，因此我们需要扩展一维DFT到二维。幸运的是，二维DFT可以通过分别在行和列上应用一维DFT来实现。

3.1 可分离性实现

二维DFT的可分离性意味着我们可以：

1. 对图像的每一行应用一维DFT
2. 对中间结果的每一列再次应用一维DFT

实现代码如下：

def dft_2d(image):
    M, N = image.shape
    # 先对每一行做一维DFT
    row_transformed = np.zeros_like(image, dtype=np.complex128)
    for i in range(M):
        row_transformed[i, :] = dft_1d(image[i, :])
    
    # 再对每一列做一维DFT
    dft_image = np.zeros_like(image, dtype=np.complex128)
    for j in range(N):
        dft_image[:, j] = dft_1d(row_transformed[:, j])
    
    return dft_image

3.2 频谱中心化

原始DFT结果的低频成分位于四角，高频在中心。为了方便观察，我们通常将低频移到中心：

def fftshift(dft_image):
    """将频域图像中心化"""
    M, N = dft_image.shape
    # 沿两个方向循环移动一半尺寸
    shifted = np.roll(dft_image, M//2, axis=0)
    shifted = np.roll(shifted, N//2, axis=1)
    return shifted

4. 图像频谱可视化与分析

现在我们可以加载真实图像并分析其频谱特性了。

4.1 加载图像并计算频谱

from skimage import data, color

# 加载示例图像并转为灰度
image = color.rgb2gray(data.astronaut())
M, N = image.shape

# 计算DFT
dft_result = dft_2d(image)
dft_shifted = fftshift(dft_result)

# 计算幅度谱（取对数便于显示）
magnitude_spectrum = np.log(np.abs(dft_shifted) + 1)  # +1避免log(0)

4.2 可视化结果

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(image, cmap='gray')
plt.title('原始图像')

plt.subplot(1, 2, 2)
plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('幅度谱')
plt.tight_layout()
plt.show()

典型频谱图会显示：

• 中心亮区域：代表图像的低频成分（整体亮度）
• 从中心向外辐射的亮线：代表图像中的边缘和纹理（高频）

4.3 与NumPy FFT对比

为了验证我们的实现，可以将其与NumPy内置的FFT进行比较：

# 使用NumPy的FFT
np_fft = np.fft.fft2(image)
np_shifted = np.fft.fftshift(np_fft)
np_magnitude = np.log(np.abs(np_shifted) + 1)

# 计算差异
difference = np.abs(magnitude_spectrum - np_magnitude)
print(f"最大差异: {np.max(difference):.6f}")

你会发现两者结果几乎相同，最大差异通常小于1e-10，验证了我们实现的正确性。

5. 实际应用与性能优化

虽然我们的实现正确，但直接使用矩阵乘法的DFT计算复杂度为O(N²)，而FFT（快速傅里叶变换）算法可以将其降低到O(N log N)。

5.1 理解FFT的优势

让我们比较两者的速度差异：

import time

# 测试不同尺寸下的计算时间
sizes = [32, 64, 128, 256]
dft_times = []
fft_times = []

for size in sizes:
    test_img = np.random.rand(size, size)
    
    # 测试DFT
    start = time.time()
    dft_2d(test_img)
    dft_times.append(time.time() - start)
    
    # 测试FFT
    start = time.time()
    np.fft.fft2(test_img)
    fft_times.append(time.time() - start)

# 绘制比较图
plt.figure(figsize=(10, 5))
plt.plot(sizes, dft_times, 'o-', label='我们的DFT实现')
plt.plot(sizes, fft_times, 's-', label='NumPy FFT')
plt.xlabel('图像尺寸')
plt.ylabel('计算时间(秒)')
plt.legend()
plt.title('DFT与FFT性能比较')
plt.grid(True)
plt.show()

随着图像尺寸增大，DFT的计算时间呈平方级增长，而FFT仅呈线性增长（在对数尺度下）。

5.2 频域滤波示例

理解傅里叶变换的真正价值在于频域滤波。例如��我们可以尝试去除图像中的高频噪声：

# 添加高频噪声
noisy = image + 0.1 * np.random.randn(*image.shape)

# 计算噪声图像的DFT
noisy_dft = dft_2d(noisy)
noisy_shifted = fftshift(noisy_dft)

# 创建低通滤波器（保留中心区域）
crow, ccol = M//2, N//2
mask = np.zeros((M, N))
r = 30  # 保留半径
mask[crow-r:crow+r, ccol-r:ccol+r] = 1

# 应用滤波器
filtered_dft = noisy_shifted * mask

# 反变换回空间域
filtered_shifted = np.fft.ifftshift(filtered_dft)
filtered_image = np.fft.ifft2(filtered_shifted).real

# 可视化结果
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(noisy, cmap='gray')
plt.title('带噪声图像')

plt.subplot(1, 2, 2)
plt.imshow(filtered_image, cmap='gray')
plt.title('频域滤波后')
plt.tight_layout()
plt.show()

你会看到高频噪声被有效去除，虽然图像变得稍微模糊（丢失了一些细节）。

6. 完整代码与扩展建议

以下是本文涉及的所有核心函数整合：

import numpy as np
import matplotlib.pyplot as plt
from skimage import data, color

def dft_1d(signal):
    N = len(signal)
    n = np.arange(N)
    k = n.reshape((N, 1))
    W = np.exp(-2j * np.pi * k * n / N)
    return np.dot(W, signal)

def dft_2d(image):
    M, N = image.shape
    row_transformed = np.zeros_like(image, dtype=np.complex128)
    for i in range(M):
        row_transformed[i, :] = dft_1d(image[i, :])
    
    dft_image = np.zeros_like(image, dtype=np.complex128)
    for j in range(N):
        dft_image[:, j] = dft_1d(row_transformed[:, j])
    
    return dft_image

def fftshift(dft_image):
    M, N = dft_image.shape
    shifted = np.roll(dft_image, M//2, axis=0)
    shifted = np.roll(shifted, N//2, axis=1)
    return shifted

# 示例使用
image = color.rgb2gray(data.astronaut())[:256, :256]  # 裁剪为256x256
dft_result = dft_2d(image)
dft_shifted = fftshift(dft_result)
magnitude_spectrum = np.log(np.abs(dft_shifted) + 1)

plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1); plt.imshow(image, cmap='gray'); plt.title('原始图像')
plt.subplot(1, 2, 2); plt.imshow(magnitude_spectrum, cmap='gray'); plt.title('幅度谱')
plt.show()

扩展建议：

1. 尝试实现快速傅里叶变换(FFT)的递归版本
2. 探索不同的频域滤波器（高通、带通、陷波）
3. 将傅里叶变换应用于音频信号处理
4. 学习短时傅里叶变换(STFT)用于时频分析

通过这次从零实现傅里叶变换的旅程，你会发现那些看似复杂的数学概念，当用代码具象化后，其实非常直观和有趣。