【问题标题】:Using FFT for 3D array representation of 2D field使用 FFT 表示 2D 场的 3D 数组
【发布时间】:2018-03-12 02:03:27
【问题描述】:

我需要得到一个复数域的傅里叶变换。我正在使用 python。

我的输入是 xy 平面中电场的 2D 快照。

我目前有一个 3D 数组 F[x][y][z] 其中 F[x][y][0] 包含实部,F[x][y]1 包含复数领域。

我当前的代码很简单,就是这样做的:

result=np.fft.fftn(F)
result=np.fft.fftshift(result)

我有以下问题:

1) 这是否正确计算了场的傅立叶变换,还是应该将场作为 2D 矩阵输入,每个元素都包含实部和虚部?

2) 我仅使用实数倍数输入了字段的复分量值(即,如果复数值是 6i,我输入了 6),这是正确的还是应该作为复数值输入(即输入为 ' 6j')?

3) 由于这在技术上是一个 2D 输入字段,我应该使用 np.fft.fft2 代替吗?这样做意味着输出不在中间。

4) 输出看起来不像我期望的 F 的傅立叶变换的样子,我不确定我做错了什么。

完整示例代码:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

x, y = np.meshgrid(np.linspace(-1,1,100), np.linspace(-1,1,100))
d = np.sqrt(x*x+y*y)
sigma, mu = .35, 0.0
g1 = np.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )

F=np.empty(shape=(300,300,2),dtype=complex)
for x in range(0,300):
        for y in range(0,300):
            if y<50 or x<100 or y>249 or x>199:
                    F[x][y][0]=g1[0][0]
                    F[x][y][1]=0j
            elif y<150:
                    F[x][y][0]=g1[x-100][y-50]
                    F[x][y][1]=0j
            else:
                    F[x][y][0]=g1[x-100][y-150]
                    F[x][y][1]=0j

F_2D=np.empty(shape=(300,300))
for x in range(0,300):
    for y in range(0,300):
            F_2D[x][y]=np.absolute(F[x][y][0])+np.absolute(F[x][y][1])

plt.imshow(F_2D)
plt.show()

result=np.fft.fftn(F)
result=np.fft.fftshift(result)

result_2D=np.empty(shape=(300,300))
for x in range(0,300):
    for y in range(0,300):
            result_2D[x][y]=np.absolute(result[x][y][0])+np.absolute(result[x][y][1])

plt.imshow(result_2D)
plt.show()

绘制 F 给出了这个:

使用 np.fft.fftn,最后显示的图像是:

使用 np.fft.fft2:

这些看起来都不像我所期望的 F 的傅立叶变换。

【问题讨论】:

  • 你应该使用F[x][y][0]而不是F[x,y,0]
  • 为此更新了我的代码,输出没有改变

标签: python numpy fft


【解决方案1】:

我在这里添加另一个答案,适合添加的代码。

答案仍然是np.fft.fft2()。这是一个例子。我稍微修改了代码。为了验证我们是否需要fft2,我丢弃了其中一个 blob,然后我们知道单个高斯 blob 应该转换为高斯 blob(具有特定相位,绘制绝对值时未显示)。我还降低了标准偏差,这样频率响应就会变宽一些。

代码:

import numpy as np
import matplotlib.pyplot as plt

x, y = np.meshgrid(np.linspace(-1,1,100), np.linspace(-1,1,100))
d = np.sqrt(x**2+y**2)
sigma, mu = .1, 0.0
g1 = np.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )
N = 300
positions = [ [150,100] ]#, [150,200] ]
sz2 = [int(x/2) for x in g1.shape]
F_2D = np.zeros([N,N])
for x0,y0 in positions:
    F_2D[ x0-sz2[0]: x0+sz2[0], y0-sz2[1]:y0+sz2[1] ] = g1 + 1j*0.

result = np.fft.fftshift(np.fft.fft2(F_2D))

plt.subplot(211); plt.imshow(F_2D)
plt.subplot(212); plt.imshow(np.absolute(result))
plt.title('$\sigma$=.1')
plt.show()

结果:

要回到原来的问题,我们只需要改变

positions = [ [150,100] , [150,200] ]sigma=.35 而不是sigma=.1

【讨论】:

  • 酷,非常感谢。只是为了仔细检查,我的错误是我应该在进行傅里叶变换之前将数组转换为二维数组?输出仍然不像我所期望的输入的傅立叶变换看起来像......:/
  • 是的,你应该给fft2一个二维(复杂)矩阵。关于您的预期结果,它是什么?也许你忘了用与序列长度相关的一些常数来缩放宽度(这是尝试计算 DTFT/连续 FT 时的常见错误,但我们只能计算 DFT)
  • 如果我将第二个 blob 的符号更改为负数(请参阅link),我希望在此链接中获得类似于 TEM1,0 图像的内容:link跨度>
  • 也许他们的距离太远了?如果我将他们的中心更改为 145 美元和 155 美元(而不是 100 美元和 200 美元),我会得到这种模式。
【解决方案2】:

您应该使用复杂的 numpy 变量(通过使用1j)并使用fft2。例如:

N = 16
x0 = np.random.randn(N,N,2)
x = x0[:,:,0] + 1j*x0[:,:,1]
X = np.fft.fft2(x)

x0 上使用fftn 将执行3D FFT,使用fft 将执行矢量方向的1D FFT。

【讨论】:

  • 这是我最初的假设——但我的类高斯函数在图像的顶部和底部(即使在应用 fftshift 之后)而不是在图像的中心有一些强度,除非我使用fftn.想法?
  • 你的空间函数是二维复高斯函数吗?在这种情况下,您应该在频域中得到一个复高斯。你能添加示例代码吗?
  • 添加了示例代码,请注意,虽然在上面的示例中 F 的虚部为零,但我打算稍后添加一个非零复数分量
猜你喜欢
  • 1970-01-01
  • 2012-07-05
  • 1970-01-01
  • 2021-12-30
  • 1970-01-01
  • 2021-09-21
  • 2018-11-13
  • 1970-01-01
相关资源
最近更新 更多