如何在Python中扩展1D FFT代码以计算图像的FFT（2D）？

Question

我正在尝试扩展一个在Python中对一维数组有效的快速傅里叶变换（FFT）代码，以便可以处理图像。实际上，我知道问题出在扩展的逻辑上。我对FFT了解不多，而我需要提交图像处理的作业。如果能给我一些提示或解决方案，我将非常感激。

这是我的代码。实际上，我正在尝试为Python创建一个FFT模块，之前在处理一维数据时已经顺利运行，得到了来自Rosetta Code网站的帮助。

from cmath import exp, pi
from math import log, ceil

def fft(f):
    N = len(f)
    if N <= 1: return f
    even = fft(f[0::2])
    odd =  fft(f[1::2])
    return [even[k] + exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)] + \
           [even[k] - exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]

def pad(f):
    n = len(f)
    N = 2 ** int(ceil(log(n, 2)))
    F = f + [0] * (N - n)
    return F, n 

def unpad(F, n):
    return F[0 : n]

def pad2(f):
    m, n = len(f), len(f[0])
    M, N = 2 ** int(ceil(log(m, 2))), 2 ** int(ceil(log(n, 2)))
    F = [ [0]*N for _ in xrange(M) ]
    for i in range(0, m):
        for j in range(0, n):
            F[i][j] = f[i][j]
    return F, m, n

def fft1D(f):
    Fu, n = pad(f)
    return fft(Fu), n

def fft2D(f):
    F, m, n = pad2(f)
    M, N = len(F), len(F[0])
    Fuv = [ [0]*N for _ in xrange(M) ]
    for i in range(0, M):
        Fxv = fft(F[i])
        for j in range(0, N):
            Fuv[i][j] = (fft(Fxv))[j]
    return Fuv, [m, n]

我用以下代码调用了这个模块：

from FFT import *
f= [0, 2, 3, 4]
F = fft1D(f)
print f, F
X, s = fft2D([[1,2,1,1],[2,1,2,2],[0,1,1,0], [0,1,1,1]])
for i in range(0, len(X)):
    print X[i]

它的输出是：

[0, 2, 3, 4] ([(9+0j), (-3+2j), (-3+0j), (-3-2j)], 4)
[(4+0j), (4-2.4492935982947064e-16j), (4+0j), (8+2.4492935982947064e-16j)]
[(8+0j), (8+2.4492935982947064e-16j), (8+0j), (4-2.4492935982947064e-16j)]
[0j, -2.33486982377251e-16j, (4+0j), (4+2.33486982377251e-16j)]
[0j, (4+0j), (4+0j), (4+0j)]

对于一维数据的结果是正确的，因为我用Matlab的输出进行了验证，但对于第二个结果，Matlab的输出是：

>> fft([1,2,1,1;2,1,2,2;0,1,1,0;0,1,1,1])

ans =

   3.0000             5.0000             5.0000             4.0000          
   1.0000 - 2.0000i   1.0000                  0 - 1.0000i   1.0000 - 1.0000i
  -1.0000             1.0000            -1.0000            -2.0000          
   1.0000 + 2.0000i   1.0000                  0 + 1.0000i   1.0000 + 1.0000i

输出结果不一样，这意味着我的代码逻辑有问题。请帮帮我，因为我到现在还没有正式学习过FFT，所以对数学部分理解得不够透彻，也许等我学完后能找到问题所在。

Answer 1

我同意isedev的看法，你应该使用numpy。它已经有一个很棒的fft（快速傅里叶变换）工具包，可以进行多维变换。

http://docs.scipy.org/doc/numpy/reference/routines.fft.html

http://docs.scipy.org/doc/numpy-1.4.x/reference/generated/numpy.fft.fft.html

Answer 2

你的代码有点难懂，但看起来你在两次操作中都是沿着同一个方向进行快速傅里叶变换（FFT）。查一下傅里叶变换的积分形式，你会发现 x 和 y 的积分是相互独立的。也就是说（抱歉，这个符号表示得不好，' 表示傅里叶空间中的一个函数）

FT(f(x, y), x) -> f'(k, y)
FT(f'(k, y), y) -> f''(k, w)

所以你需要做的是对每一行进行快速傅里叶变换（也就是做 N 次一维的 FFT），然后把结果放到一个新的数组里（这一步是把 f(x, y) -> f'(k, y)）。接着，再对这个结果数组的每一列进行快速傅里叶变换（也就是做 M 次一维的 FFT），然后把这些结果放到另一个新的数组里（这一步是把 f'(k, y) -> f''(k, w)）。