别怕数学用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_domain2.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, d1/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的可分离性意味着我们可以对图像的每一行应用一维DFT对中间结果的每一列再次应用一维DFT实现代码如下def dft_2d(image): M, N image.shape # 先对每一行做一维DFT row_transformed np.zeros_like(image, dtypenp.complex128) for i in range(M): row_transformed[i, :] dft_1d(image[i, :]) # 再对每一列做一维DFT dft_image np.zeros_like(image, dtypenp.complex128) for j in range(N): dft_image[:, j] dft_1d(row_transformed[:, j]) return dft_image3.2 频谱中心化原始DFT结果的低频成分位于四角高频在中心。为了方便观察我们通常将低频移到中心def fftshift(dft_image): 将频域图像中心化 M, N dft_image.shape # 沿两个方向循环移动一半尺寸 shifted np.roll(dft_image, M//2, axis0) shifted np.roll(shifted, N//2, axis1) return shifted4. 图像频谱可视化与分析现在我们可以加载真实图像并分析其频谱特性了。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, cmapgray) plt.title(原始图像) plt.subplot(1, 2, 2) plt.imshow(magnitude_spectrum, cmapgray) 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-, labelNumPy 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:crowr, ccol-r:ccolr] 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, cmapgray) plt.title(带噪声图像) plt.subplot(1, 2, 2) plt.imshow(filtered_image, cmapgray) 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, dtypenp.complex128) for i in range(M): row_transformed[i, :] dft_1d(image[i, :]) dft_image np.zeros_like(image, dtypenp.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, axis0) shifted np.roll(shifted, N//2, axis1) 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, cmapgray); plt.title(原始图像) plt.subplot(1, 2, 2); plt.imshow(magnitude_spectrum, cmapgray); plt.title(幅度谱) plt.show()扩展建议尝试实现快速傅里叶变换(FFT)的递归版本探索不同的频域滤波器高通、带通、陷波将傅里叶变换应用于音频信号处理学习短时傅里叶变换(STFT)用于时频分析通过这次从零实现傅里叶变换的旅程你会发现那些看似复杂的数学概念当用代码具象化后其实非常直观和有趣。
别怕数学！用Python从零实现图像傅里叶变换（附完整代码与频谱图分析）
别怕数学用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_domain2.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, d1/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的可分离性意味着我们可以对图像的每一行应用一维DFT对中间结果的每一列再次应用一维DFT实现代码如下def dft_2d(image): M, N image.shape # 先对每一行做一维DFT row_transformed np.zeros_like(image, dtypenp.complex128) for i in range(M): row_transformed[i, :] dft_1d(image[i, :]) # 再对每一列做一维DFT dft_image np.zeros_like(image, dtypenp.complex128) for j in range(N): dft_image[:, j] dft_1d(row_transformed[:, j]) return dft_image3.2 频谱中心化原始DFT结果的低频成分位于四角高频在中心。为了方便观察我们通常将低频移到中心def fftshift(dft_image): 将频域图像中心化 M, N dft_image.shape # 沿两个方向循环移动一半尺寸 shifted np.roll(dft_image, M//2, axis0) shifted np.roll(shifted, N//2, axis1) return shifted4. 图像频谱可视化与分析现在我们可以加载真实图像并分析其频谱特性了。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, cmapgray) plt.title(原始图像) plt.subplot(1, 2, 2) plt.imshow(magnitude_spectrum, cmapgray) 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-, labelNumPy 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:crowr, ccol-r:ccolr] 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, cmapgray) plt.title(带噪声图像) plt.subplot(1, 2, 2) plt.imshow(filtered_image, cmapgray) 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, dtypenp.complex128) for i in range(M): row_transformed[i, :] dft_1d(image[i, :]) dft_image np.zeros_like(image, dtypenp.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, axis0) shifted np.roll(shifted, N//2, axis1) 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, cmapgray); plt.title(原始图像) plt.subplot(1, 2, 2); plt.imshow(magnitude_spectrum, cmapgray); plt.title(幅度谱) plt.show()扩展建议尝试实现快速傅里叶变换(FFT)的递归版本探索不同的频域滤波器高通、带通、陷波将傅里叶变换应用于音频信号处理学习短时傅里叶变换(STFT)用于时频分析通过这次从零实现傅里叶变换的旅程你会发现那些看似复杂的数学概念当用代码具象化后其实非常直观和有趣。
