Python 实现二维 DFT 时结果异常，如何修正？

2026-11-04数据阿明

结论

主要错误是把原图各行的一维 DFT 结果与各列的一维 DFT 结果相加，当作二维 DFT。正确做法是先对每一行做完整的一维 DFT，再对得到的复数矩阵逐列做一维 DFT，将两次变换串接起来。

代码中的指数项还缺少虚数单位，也没有包含完整的频率索引与采样位置的乘积，因此计算结果不是 DFT。

给出的红色通道也并非全为 255：有 7 个像素的 R 值为 250。即使修正算法，这份数据的红色频谱仍会包含非零的非直流分量。

二维 DFT 应该怎样计算？

设图像单个通道为 f[y][x]，高度为 H，宽度为 W。采用正变换不归一化的约定：

[ F[u,v]=\sum_{y=0}^{H-1}\sum_{x=0}^{W-1} f[y,x]\exp\left(-2\pi j\left(\frac{uy}{H}+\frac{vx}{W}\right)\right). ]

其中，x、y 是空间位置，u、v 是频率索引，j 是虚数单位。

利用指数函数的乘法关系，可以先计算：

[ R[y,v]=\sum_{x=0}^{W-1}f[y,x]\exp\left(-2\pi j\frac{vx}{W}\right), ]

再计算：

[ F[u,v]=\sum_{y=0}^{H-1}R[y,v]\exp\left(-2\pi j\frac{uy}{H}\right). ]

第一步为每行产生 W 个频率系数，结果保留为一个 H×W 复数矩阵；第二步变换这个矩阵的每一列。先列后行也可以。

因此，“一行得到 1、一列得到 2，所以交点为 3”这个理解不成立。每个二维频率系数都要结合所有像素计算。

原代码有哪些具体错误？

指数项没有虚数单位。

原代码使用：
```
cmath.exp(-1*cmath.pi*i/8)
```
这里的 -1 是实数，不是虚数单位。实数指数产生实数权重，即使通过 cmath 返回复数类型，虚部仍然为零。这解释了当前输出中大量的 +0j。

一维 DFT 对应的 Python 表达式应为：
```
cmath.exp(-2j * cmath.pi * k * n / N)
```
k 是输出频率索引，n 是输入采样位置。原代码也没有完整包含 k*n 这个乘积。
每行只保存了一个数。

长度为 16 的一维 DFT 应产生 16 个复数系数。fourierRowsRed[i] 只保存一个标量，无法表示整行频谱。行变换的中间结果也应是二维矩阵。
所谓列循环仍在访问行。

固定 i、改变 j 时，redArray[i][j] 遍历的是第 i 行。若要遍历原图的一列，应固定列索引、改变行索引；在可分离实现的第二步，则应遍历行变换结果的一列。
相加不能完成第二次变换。
```
redArrayFourier[i][j] = fourierRowsRed[i]+fourierColumnsRed[j]
```
这一步没有二维 DFT 的数学依据。应将行变换结果作为列变换的输入。
延迟写入会漏掉最后一次累加。

按描述推断的循环缩进，if i>0 在下一轮才保存上一轮的累加值。最后一轮结束后没有后续写入，数组最后一项因此保留初始值。可以在每个频率系数计算完成后直接赋值，避免这个问题。

原始粘贴代码还缺少 Python 必需的缩进，运行文件时需要保留正确缩进。以上分析针对描述和输出所体现的计算逻辑。

如何替换 DFT 计算部分？

保留文件解析得到的 redArray、greenArray 和 blueArray，用下面的代码替换原来的行列累加与相加部分。这里假设每个数组都是非空、各行等长的二维矩阵。

import cmath


def dft1d(values):
    N = len(values)
    return [
        sum(
            values[n] * cmath.exp(-2j * cmath.pi * k * n / N)
            for n in range(N)
        )
        for k in range(N)
    ]


def dft2d(image):
    H = len(image)
    W = len(image[0])

    # First transform each row, retaining every frequency coefficient.
    row_spectrum = [dft1d(row) for row in image]

    # Then transform each column of the intermediate complex matrix.
    spectrum = [[0j for _ in range(W)] for _ in range(H)]
    for v in range(W):
        column = [row_spectrum[y][v] for y in range(H)]
        transformed_column = dft1d(column)
        for u in range(H):
            spectrum[u][v] = transformed_column[u]

    return spectrum


redArrayFourier = dft2d(redArray)
greenArrayFourier = dft2d(greenArray)
blueArrayFourier = dft2d(blueArray)

print(redArrayFourier)
print(greenArrayFourier)
print(blueArrayFourier)

这份实现直接使用 DFT 定义，便于检查频率和位置索引。RGB 三个通道分别计算；第二次变换要保留第一次产生的复数值，不要提前取实部或绝对值。

怎样判断红色和绿色通道的结果？

采用正变换不归一化的约定时，直流分量就是所有像素值之和：

[ F[0,0]=\sum_y\sum_x f[y,x]. ]

如果红色通道确实全部为 255，则：

[ F_R[0,0]=16\times16\times255=65280, ]

其他频率系数理论上为零。

但提供的数据中，第 4 行有 5 个 R 值为 250 的像素，第 13 行有 2 个，这里按自然行号从 1 开始计数。因此，这份数据对应的红色直流分量应为：

[ F_R[0,0]=65280-7\times5=65245. ]

这些偏离 255 的像素也会产生非直流频率分量。可以先检查实际解码出的数组：

print("Red values:", sorted({value for row in redArray for value in row}))
print("Expected DC:", sum(sum(row) for row in redArray))
print("Computed DC:", redArrayFourier[0][0])

若文件内容与题目给出的数据一致，红色取值应包含 250 和 255，像素总和应为 65245。这是根据输入数据算出的预期值，不是实际运行结果。

绿色通道全部为零，其 DFT 也应全部为零。不过，错误算法同样可能把全零输入映射为全零，所以不能只凭绿色通道的结果判断实现是否正确。

还可以单独构造常量输入来检查算法：

constant_red = [[255 for _ in range(16)] for _ in range(16)]
F = dft2d(constant_red)

tolerance = 1e-6
assert abs(F[0][0] - 65280) < tolerance
assert all(
    abs(F[u][v]) < tolerance
    for u in range(16)
    for v in range(16)
    if (u, v) != (0, 0)
)

浮点计算可能留下很小的残差，因此应使用容差判断零值，不宜要求所有非直流系数都严格等于 0j。

如何查看图像频谱？

DFT 输出的是复数系数，不能直接当作普通 RGB 像素显示。查看频谱时，通常先取幅度，再用对数压缩显示范围，并把零频移动到中心。

如果安装了 NumPy 和 Matplotlib，可以用下面的代码显示上面计算出的蓝色通道频谱：

import numpy as np
import matplotlib.pyplot as plt

F_blue = np.asarray(blueArrayFourier, dtype=complex)
display_spectrum = np.log1p(np.abs(np.fft.fftshift(F_blue)))

plt.imshow(display_spectrum, cmap="gray")
plt.title("Blue channel spectrum")
plt.colorbar()
plt.show()

fftshift 只调整显示位置，不会修正变换结果。原始频谱中的直流位置是 [0][0]，移位后才位于图像中心。幅度图适合观察频率强弱；若要重建图像，还需要保留完整的复数系数及其相位。