3
我试图创建一个高斯随机场,通过在傅里叶空间中创建一个网格,然后通过反傅里叶变换来获得随机场。为此,逆傅立叶变换图像需要是实值的。我似乎在10^-18 - -22阶的网格的虚部中得到了残差,所以我期望这是FFT中的数值误差。尽管如此,图像的实际部分在像素上显示了一个奇怪的棋盘图案,其中像素从正面跳到负面。为了查看FFT功能是否正常工作,我尝试转换一个高斯,这应该返回另一个高斯,并且图像中存在棋盘图案。当获取图像的绝对值时,它看起来很好,但我也需要它来允许我的高斯随机场的负值。FFT后的棋盘图案
对于高斯我使用下面的代码的傅立叶变换:
#! /usr/bin/env python
import numpy as n
import math as m
import pyfits
def fourierplane(a):
deltakx = 2*a.kxmax/a.dimkx #stepsize in k_x
deltaky = 2*a.kymax/a.dimky #stepsize in k_y
plane = n.zeros([a.dimkx,a.dimky]) #empty matrix to be filled in for the Fourier grid
for y in range(n.shape(plane)[0]):
for x in range(n.shape(plane)[1]):
#Defining coordinates centred at x = N/2, y = N/2
i1 = x - a.dimkx/2
j1 = y - a.dimky/2
#creating values to fill in in the grid:
kx = deltakx*i1 #determining value of k_x at gridpoint
ky = deltaky*j1 #determining value of k_y at gridpoint
k = m.sqrt(kx**2 + ky**2) #magnitude of k-vector
plane[y][x] = m.e**(-(k**2)/(2*a.sigma_k**2)) #gaussian
return plane
def substruct():
class fougrid:
pass
grid = fougrid()
grid.kxmax = 2.00 #maximum value k_x
grid.kymax = 2.00 #maximum value k_y
grid.sigma_k = (1./20.)*grid.kxmax #width of gaussian
grid.dimkx = 1024
grid.dimky= 1024
fplane = fourierplane(grid) #creating the Fourier grid
implane = n.fft.ifftshift(n.fft.ifft2(fplane)) #inverse Fourier transformation of the grid to get final image
##################################################################
#seperating real and imaginary part of the Fourier transformed grid
##################################################################
realimplane = implane.real
imagimplane = implane.imag
#taking the absolute value:
absimplane = n.zeros(n.shape(implane))
for a in range(n.shape(implane)[0]):
for b in range(n.shape(implane)[1]):
absimplane[a][b] = m.sqrt(implane[a][b].real**2 + implane[a][b].imag**2)
#saving images to files:
pyfits.writeto('randomfield.fits',realimplane) #real part of the image grid
pyfits.writeto('fplane.fits',fplane) #grid in fourier space
pyfits.writeto('imranfield.fits',imagimplane) #imaginary part of the image grid
pyfits.writeto('absranfield.fits',absimplane) #real part of the image grid
substruct() #running the script
有没有人有任何想法,这种模式是如何创建的,以及如何解决这个问题?
解决了它。 在代码中,上fplane进行ifft2之前,fplane需要被移位为好,这样: implane = n.fft.ifftshift(n.fft.ifft2(n.fft.fftshift(fplane)) ) – Mizuti 2012-04-19 12:16:53