当前位置： 首页 > news >正文

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

news 2026/5/31 15:21:55

别怕数学！用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)用于时频分析

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

http://www.cnnetsun.cn/news/2521122.html

相关文章：

  • 告别训练慢和显存焦虑：RTMDet实战中那些你没注意到的工程优化细节（附代码）
  • AXI总线安全访问机制与寄存器布局实践
  • C语言高级笔记
  • Keil C51递归调用警告处理与工程配置详解
  • ARM嵌入式开发中DS-5内存优化与JVM调优实战
  • 大麦网自动化抢票解决方案：告别手动抢票的低效困境
  • fuckZHS：智慧树课程自动化学习脚本深度解析与逆向工程技术实现
  • 可以快速引蜘蛛的蜘蛛池是什么？
  • Webdash API详解：如何通过RESTful接口扩展和集成外部系统
  • Zhui组件库开发指南：从环境搭建到贡献代码的完整路线图
  • Beat Saber版本管理终极解决方案：BSManager完全指南
  • 3分钟搞定系统镜像烧录！Balena Etcher：开源免费的跨平台烧录神器
  • Ventoy主题定制完全指南：让你的启动界面焕然一新！
  • Scribd电子书离线下载：构建个人数字图书馆的一站式自动化解决方案
  • “冠珠·美乐童行”公益行动走进广州市增城区高滩小学，唱响爱、筑就美
  • sdk-manager-plugin历史与演进：从诞生到废弃的完整技术演进路线图
  • 3个真实场景揭秘：res-downloader如何帮你节省90%的视频收集时间
  • 城市交通气候适应：从生物滞留池到透水铺装的工程实践
  • 3D高斯泼溅技术实现实时4D天气模拟
  • 均衡传播算法(EP)原理与硬件实现优势
  • 微信小程序 零工市场服务系统
  • 量子退火与组合优化：LDA框架的创新应用
  • Linux服务与权限安全加固——从“服务起不来“到“安全合规“的5层防御体系
  • 《Sysinternals实战指南》ZoomIt 学习笔记（11.10）：键入模式——在桌面上直接打字讲解的最佳实践
  • 为什么选择SecHex-Spoofy？对比5款HWID工具，这款开源神器究竟强在哪里
  • Recipe协议：基于TEE的BFT复制协议设计与优化
  • AI INFRA之NVIDIA GPUDirect节点内和节点间通信原理详解
  • 计算机视觉——九、图像分割
  • PHP 的 resource（如数据库连接、文件句柄）不能被序列化。
  • H3CSE 高性能园区网：生成树保护机制