如何将自行实现的1D DFT/FFT算法应用于灰度图像的傅里叶变换

核心原理 #

2D离散傅里叶变换可以通过两次1D变换实现：先对图像的每一行做1D傅里叶变换，再对得到的结果的每一列做1D傅里叶变换，即可得到完整的2D变换结果，不需要重新修改你已有的1D DFT/FFT函数。

前置注意事项 #

• 你当前实现的基2 FFT仅支持输入长度为2的整数次幂，处理图像前需要将图像的宽、高调整为2的整数次幂，可通过裁剪或边缘补零实现，否则FFT函数会报错
• 输入图像需要转换为浮点型计算，避免uint8类型的精度丢失
• 可视化频谱前需要执行fftshift操作将低频分量移到图像中心，同时对幅度谱做对数拉伸解决直流分量亮度过高、其他分量无法显示的问题

新增2D变换代码 #

你可以直接在现有代码基础上添加以下两个函数：

def DFT2D(img):
    rows, cols = img.shape
    # 先处理每一行
    row_res = ny.zeros_like(img, dtype=complex)
    for i in range(rows):
        row_res[i, :] = DFT(img[i, :])
    # 再处理每一列
    col_res = ny.zeros_like(row_res, dtype=complex)
    for j in range(cols):
        col_res[:, j] = DFT(row_res[:, j])
    return col_res

def FFT2D(img):
    rows, cols = img.shape
    # 先处理每一行
    row_res = ny.zeros_like(img, dtype=complex)
    for i in range(rows):
        row_res[i, :] = FFT(img[i, :])
    # 再处理每一列
    col_res = ny.zeros_like(row_res, dtype=complex)
    for j in range(cols):
        col_res[:, j] = FFT(row_res[:, j])
    return col_res

# 补充fftshift实现，用于调整频谱显示
def fftshift2d(img):
    rows, cols = img.shape
    return ny.roll(ny.roll(img, rows//2, axis=0), cols//2, axis=1)

灰度图像处理示例代码 #

# 读入灰度图像，替换为你自己的图像路径
img = mplot.imread('gray_image.jpg')
# 若读入为三通道彩色图，可转灰度：img = ny.mean(img, axis=2)
# 调整图像尺寸为2的整数次幂，示例取512*512
img = ny.pad(img, ((0, 512-img.shape[0]), (0, 512-img.shape[1])), mode='constant')
img = img.astype(ny.float64)

# 计算2D FFT
fft2d_res = FFT2D(img)
# 计算幅度谱和相位谱
mag_spectrum = ny.abs(fft2d_res)
phase_spectrum = ny.angle(fft2d_res)

# 可视化
mplot.figure(figsize=(12,6))
mplot.subplot(121)
# 幅度谱做对数拉伸+中心平移
mplot.imshow(fftshift2d(ny.log(1 + mag_spectrum)), cmap='gray')
mplot.title('FFT幅度谱')
mplot.axis('off')

mplot.subplot(122)
mplot.imshow(fftshift2d(phase_spectrum), cmap='gray')
mplot.title('FFT相位谱')
mplot.axis('off')
mplot.show()

其他优化提示 #

• 你原代码中相位谱的计算用Y/abs(Y)会出现除以零的问题，直接调用ny.angle(Y)即可得到正确的相位值
• 大尺寸图像计算DFT会非常慢，建议只用FFT做2D变换
• 如果需要实现逆变换，只需要将1D变换中的指数项改为正号，最后结果除以序列长度即可，2D逆变换同样用两次1D逆变换实现。

内容的提问来源于stack exchange，提问作者J.Smith