python中的DFT矩阵答案

【问题标题】：DFT matrix in pythonpython中的DFT矩阵
【发布时间】：2013-11-13 09:44:33
【问题描述】：

在 python 中为二维 DFT 获取 DFT matrix 的最简单方法是什么？我在numpy.fft 中找不到这样的功能。谢谢！

【问题讨论】：

【解决方案1】：

最简单也是最快的方法是使用 SciPy 中的 fft。

import scipy as sp

def dftmtx(N):
    return sp.fft(sp.eye(N))

如果您知道更快的方法（可能更复杂），我将不胜感激您的意见。

只是为了让它与主要问题更相关——你也可以用 numpy 来做：

import numpy as np

dftmtx = np.fft.fft(np.eye(N))

当我对它们进行基准测试时，我的印象是 scipy 的速度稍微快了一点，但我还没有彻底完成，而且是很久以前的事了，所以不要相信我的话。

以下是关于 python 中的 FFT 实现的非常好的来源： http://nbviewer.ipython.org/url/jakevdp.github.io/downloads/notebooks/UnderstandingTheFFT.ipynb 而是从速度的角度来看，但在这种情况下，我们实际上可以看到有时它也很简单。

【讨论】：

我问的时候它不在那里，你看。 :) 我 +1d 你，所以它会冒泡。
好主意，谢谢。很久以前，我一直在为 DFT 矩阵生成而苦苦挣扎，并花了一段时间来缩小范围以找到一个合理的解决方案。当您生成大问题时，DFT 矩阵生成是一个巨大的性能瓶颈。我会投票赞成这个问题，但我还是新手。

【解决方案2】：

我不认为这是内置的。但是，直接计算很简单：

import numpy as np
def DFT_matrix(N):
    i, j = np.meshgrid(np.arange(N), np.arange(N))
    omega = np.exp( - 2 * pi * 1J / N )
    W = np.power( omega, i * j ) / sqrt(N)
    return W

编辑对于二维 FFT 矩阵，您可以使用以下内容：

x = np.zeros(N, N) # x is any input data with those dimensions
W = DFT_matrix(N)
dft_of_x = W.dot(x).dot(W)

【讨论】：

好的，这确实是 1d DFT 的矩阵；我需要将其更改为用于 2d DFT 的 N^2xN^2 矩阵。
@UriCohen：我编辑了答案，有点描述了如何做你想做的事。如果你想要一个矩阵，你必须做一些代数:)
如何进行逆 DFT？ @AlexI
@AlexI 虽然它可能有类似的效果，但从技术上讲你应该应用它：dft_of_x = W.dot(x).dot(W.T)。你无法逃脱dft_of_x = (W.dot(W)).dot(x)。从左到右只应用一个矩阵只能变换x的列；这将是（不是 DFT）一维变换。
只有一句话.. 当我将您的 DFT_matrix 与 np.fft.fft(np.eye(N)) 中的 DFT_matrix 进行比较时，我看到 np.pi “缺失”的一个因素。我知道这只是一个比例因子，但想提一下

【解决方案3】：

截至 scipy 0.14 有一个内置的scipy.linalg.dft：

以 16 点 DFT 矩阵为例：

>>> import scipy.linalg
>>> import numpy as np
>>> m = scipy.linalg.dft(16)

验证单一属性，注释矩阵未缩放因此16*np.eye(16)：

>>> np.allclose(np.abs(np.dot( m.conj().T, m )), 16*np.eye(16))
True

对于 2D DFT 矩阵，这只是张量积的问题，或者在这种情况下特别是 Kronecker 积，因为我们正在处理矩阵代数。

>>> m2 = np.kron(m, m) # 256x256 matrix, flattened from (16,16,16,16) tensor

现在我们可以给它一个平铺的可视化，它是通过将每一行重新排列成一个方块来完成的

>>> import matplotlib.pyplot as plt
>>> m2tiled = m2.reshape((16,)*4).transpose(0,2,1,3).reshape((256,256))
>>> plt.subplot(121)
>>> plt.imshow(np.real(m2tiled), cmap='gray', interpolation='nearest')
>>> plt.subplot(122)
>>> plt.imshow(np.imag(m2tiled), cmap='gray', interpolation='nearest')
>>> plt.show()

结果（实部和虚部分别）：

如您所见，它们是二维 DFT 基函数

Link 到文档

【讨论】：

【解决方案4】：

@亚历克斯|基本上是正确的，我在这里添加我用于2-d DFT的版本：

def DFT_matrix_2d(N):
    i, j = np.meshgrid(np.arange(N), np.arange(N))
    A=np.multiply.outer(i.flatten(), i.flatten())
    B=np.multiply.outer(j.flatten(), j.flatten())
    omega = np.exp(-2*np.pi*1J/N)
    W = np.power(omega, A+B)/N
    return W

【讨论】：

【解决方案5】：

Lambda 函数也可以：

dftmtx = lambda N: np.fft.fft(np.eye(N))

您可以使用 dftmtx(N) 调用它。示例：

In [62]: dftmtx(2)
Out[62]: 
array([[ 1.+0.j,  1.+0.j],
       [ 1.+0.j, -1.+0.j]])

【讨论】：

【解决方案6】：

如果您希望将 2D DFT 计算为单个矩阵运算，则有必要将您希望在其上计算 DFT 的矩阵 X 分解为一个向量，因为 DFT 的每个输出对输入，而单个方阵乘法不具备这种能力。注意确保我们正确处理索引，我发现以下工作：

M = 16
N = 16
X = np.random.random((M,N)) + 1j*np.random.random((M,N))
Y = np.fft.fft2(X)
W = np.zeros((M*N,M*N),dtype=np.complex)

hold = []
for m in range(M):
    for n in range(N):
        hold.append((m,n))

for j in range(M*N):
    for i in range(M*N):
        k,l = hold[j]
        m,n = hold[i]
        W[j,i] = np.exp(-2*np.pi*1j*(m*k/M + n*l/N))

np.allclose(np.dot(W,X.ravel()),Y.ravel())
True

如果您希望将归一化更改为正交，您可以除以 1/sqrt(MN)，或者如果您希望进行逆变换，只需更改指数中的符号即可。

【讨论】：

【解决方案7】：

这可能有点晚了，但是使用 NumPy 的 vander 创建 DFT 矩阵有一个更好的替代方法，它执行得更快

此外，此实现不使用循环（显式）

def dft_matrix(signal):

    N = signal.shape[0]  # num of samples
    w = np.exp((-2 * np.pi * 1j) / N)  # remove the '-' for inverse fourier
    r = np.arange(N)
    w_matrix = np.vander(w ** r, increasing=True)  # faster than meshgrid
    return w_matrix

如果我没记错的话，主要的改进是这种方法从（已经计算的）之前的元素中生成了幂的元素

您可以在文档中了解vander： numpy.vander

【讨论】：