2012-04-19 83 views
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 

有没有人有任何想法,这种模式是如何创建的,以及如何解决这个问题?

+0

解决了它。 在代码中,上fplane进行ifft2之前,fplane需要被移位为好,这样: implane = n.fft.ifftshift(n.fft.ifft2(n.fft.fftshift(fplane)) ) – Mizuti 2012-04-19 12:16:53

回答

1

每当在一个DFT域中看到意外的交替符号时,可能意味着其他DFT域中的数据在阵列中间旋转(类似于fftshift)。如果在一个域中有实数值的对称“驼峰”,那么在数组元素0(而不是数组元素n/2)上对中驼峰将是最有可能不会在变换域中产生交替符号的布局。