在 python中获取二维 DFT 的DFT 矩阵的最简单方法是什么?我在numpy.fft中找不到这样的功能。谢谢!
7 回答
最简单且最可能最快的方法是使用 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 这是从速度的角度来看,但在这种情况下,我们实际上可以看到有时它也很简单。
我不认为这是内置的。但是,直接计算很简单:
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)
截至 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()
结果(实部和虚部分别):
如您所见,它们是 2D DFT 基函数
链接到文档
@亚历克斯| 基本上是正确的,我在这里添加我用于二维 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
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]])
如果您希望将 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),或者如果您希望进行逆变换,只需更改指数中的符号即可。
这可能有点晚了,但是创建 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