在NumPy中沿3D阵列的每个轴应用DFT矩阵？

Question

我可以先获得给定大小的DFT矩阵,表示n通过

import numpy as np
n = 64
D = np.fft.fft(np.eye(n))

FFT当然只是应用于D矢量的快速算法:

x = np.random.randn(n)
ft1 = np.dot(D,x)
print( np.abs(ft1 - fft.fft(x)).max() )
# prints near double precision roundoff

可以通过应用D矩阵的行和列来获得2D FFT :

x = np.random.randn(n,n)
ft2 = np.dot(x, D.T) # Apply D to rows.
ft2 = np.dot(D, ft2) # Apply D to cols.
print( np.abs(ft2 - fft.fft2(x)).max() )
# near machine round off again

我如何类似地计算三维离散傅立叶变换？

也就是说,

x = np.random.randn(n,n,n)
ft3 = # dot operations using D and x
print( np.abs(ft3 - fft.fftn(x)).max() )
# prints near zero

基本上,我认为我需要应用于卷D中的每个列向量,然后应用卷中的每个行向量,最后是每个"深度向量".但我不知道如何使用它dot.

Answer 1

您可以使用einsum表达式对每个索引执行转换:

x = np.random.randn(n, n, n)
ft3 = np.einsum('ijk,im->mjk', x, D)
ft3 = np.einsum('ijk,jm->imk', ft3, D)
ft3 = np.einsum('ijk,km->ijm', ft3, D)
print(np.abs(ft3 - np.fft.fftn(x)).max())
1.25571216554e-12

这也可以写成一个NumPy步骤:

ft3 = np.einsum('ijk,im,jn,kl->mnl', ft3, D, D, D, optimize=True)

如果没有optimize参数(在NumPy 1.12+中可用),它将会非常缓慢.您也可以使用每个步骤dot,但需要进行一些整形和转置.在NumPy 1.14+中,该einsum功能将自动检测BLAS操作并为您执行此操作.