20

在 python中获取二维 DFT 的DFT 矩阵的最简单方法是什么?我在numpy.fft中找不到这样的功能。谢谢!

4

7 回答 7

17

最简单且最可能最快的方法是使用 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 这是从速度的角度来看,但在这种情况下,我们实际上可以看到有时它也很简单。

于 2014-10-29T18:04:43.140 回答
14

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

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

编辑对于 2D 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)
于 2013-11-02T06:38:13.180 回答
7

截至 scipy0.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 基础

如您所见,它们是 2D DFT 基函数

链接到文档

于 2017-02-17T05:25:26.730 回答
2

@亚历克斯| 基本上是正确的,我在这里添加我用于二维 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
于 2013-11-02T22:05:15.780 回答
2

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]])
于 2015-06-24T00:00:46.707 回答
2

如果您希望将 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),或者如果您希望进行逆变换,只需更改指数中的符号即可。

于 2020-02-19T00:03:02.597 回答
1

这可能有点晚了,但是创建 DFT 矩阵有一个更好的选择,它使用 NumPy 的执行速度更快vander

此外,此实现不使用循环(明确)

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

于 2020-11-11T12:50:46.860 回答