将 FFT 用于 2D 场的 3D 阵列表示
Using FFT for 3D array representation of 2D field
我需要获取复杂场的傅里叶变换。我正在使用 python.
我的输入是 xy 平面中电场的二维快照。
我目前有一个 3D 数组 F[x][y][z],其中 F[x][y][0] 包含实部,F[x][y]1 包含字段的复杂组件。
我当前的代码非常简单,就是这样做的:
result=np.fft.fftn(F)
result=np.fft.fftshift(result)
我有以下问题:
1) 这是否正确计算场的傅里叶变换,或者场是否应该作为二维矩阵输入,每个元素都包含实部和虚部?
2) 我仅使用实数倍数输入了字段的复分量值(即如果复数值是 6i 我输入了 6),这是正确的还是应该作为复数值输入(即输入作为'6j')?
3) 由于这在技术上是一个二维输入字段,我应该改用 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 的傅里叶变换。
您应该使用复杂的 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
将执行 vector-wise 1D FFT。
我在这里添加另一个答案,适合添加的代码。
答案仍然是np.fft.fft2()
。这是一个例子。我稍微修改了代码。为了验证我们是否需要 fft2
我丢弃了其中一个斑点,然后我们知道单个高斯斑点应该转换为高斯斑点(具有特定相位,绘制绝对值时未显示)。我还降低了标准偏差,这样频率响应就会变宽一点。
代码:
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
.
我需要获取复杂场的傅里叶变换。我正在使用 python.
我的输入是 xy 平面中电场的二维快照。
我目前有一个 3D 数组 F[x][y][z],其中 F[x][y][0] 包含实部,F[x][y]1 包含字段的复杂组件。
我当前的代码非常简单,就是这样做的:
result=np.fft.fftn(F)
result=np.fft.fftshift(result)
我有以下问题:
1) 这是否正确计算场的傅里叶变换,或者场是否应该作为二维矩阵输入,每个元素都包含实部和虚部?
2) 我仅使用实数倍数输入了字段的复分量值(即如果复数值是 6i 我输入了 6),这是正确的还是应该作为复数值输入(即输入作为'6j')?
3) 由于这在技术上是一个二维输入字段,我应该改用 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 的傅里叶变换。
您应该使用复杂的 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
将执行 vector-wise 1D FFT。
我在这里添加另一个答案,适合添加的代码。
答案仍然是np.fft.fft2()
。这是一个例子。我稍微修改了代码。为了验证我们是否需要 fft2
我丢弃了其中一个斑点,然后我们知道单个高斯斑点应该转换为高斯斑点(具有特定相位,绘制绝对值时未显示)。我还降低了标准偏差,这样频率响应就会变宽一点。
代码:
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
.