2012-04-19 4 views
3

Я пытаюсь создать случайное поле Гаусса, создав сетку в пространстве Фурье, а затем обратный Фурье, переводя ее, чтобы получить случайное поле. Для этого обратное преобразованное Фурье изображение должно быть вещественным. Кажется, я получаю остатки в мнимой части сетки порядка 10^-18 - -22, поэтому я ожидал, что это будут числовые ошибки в БПФ. Реальная часть изображения отображает странный шаблон шахматной доски на пикселях, хотя пиксели переходят с положительного на отрицательный. Чтобы проверить, правильно ли функционирует 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

Решил. В коде, до того, как ifft2 выполняется на fplane, fplane должна быть смещена, так что: implane = n.fft.ifftshift (n.fft.ifft2 (n.fft.fftshift (fplane))) – Mizuti

ответ

1

Всякий раз, когда вы видите неожиданные знакопеременные знаки в одном домене DFT, это может означать, что данные в другой области DFT были повернуты наполовину через массив (аналогично fftshift). Если у вас есть симметричный «горб» реальных значений в одном домене, то центрирование этого горба на элементе массива 0 (вместо элемента массива n/2) будет состоянием, которое, скорее всего, не приведет к появлению чередующихся знаков в домене преобразования.

Смежные вопросы