Ось частоты в Numpy fft

Question

У меня есть изображение размером 32x32 пикселя, которое я хотел бы вычислить 1D спектр мощности (усредненный вдоль оси y изображения). Вот что я делаю.

import numpy as np 

size_patch=32 

# Take the fourier transform of the image. 
F1_obs = np.fft.fft2(image_obs) 

# Now shift the quadrants around so that low spatial frequencies are in 
# the center of the 2D fourier transformed image. 
F2_obs = np.fft.fftshift(F1_obs) 

# Calculate a 2D power spectrum 
psd2D_obs=np.abs(F2_obs)**2 

freq = np.fft.fftfreq(size_patch, d=1.) 

#Compute the 1D power spectrum (averaged along y axis) 
psd1D_obs[i] = np.average(psd2D_obs,axis = 0)[size_patch/2:] # we keep the end values of the array : the fft is symmetric 

У меня есть некоторые проблемы с пониманием того, что точно построено как ось х в спектре мощности. Является ли это волновым числом или пространственной частотой? Какова конвенция, принятая здесь? Является ли естественная единица циклов/пикселей? Документ для numpy.fft.fftfreq слишком расплывчатый для меня. Это может быть очень простой вопрос, но я ничего не нашел в любом месте. Не могли бы вы просветить меня?

Answer 1

Если вы хотите усреднить 2D-спектр, вы должны сделать радиальное среднее значение, а не вдоль оси. Радиальное среднее между r и r + dr представляет собой сумму значений пикселей, которые расстояние r до центра находится между r и dr, деленное на количество этих пикселей. Это возможно с помощью цикла, но я бы предпочел использовать numpy.histogram дважды.

Вот класс, который выполняет эту работу, с возможностью оконной обработки. Он оптимизирован для обычного случая, когда у вас много изображений (или патчей) одного размера для лечения. Я также использую numexpr, чтобы ускорить то, что может быть.

В дальнейшем я называю свою пространственную частоту q.

import numpy as np 
import numexpr 

class ImageStructureFactor: 
    """A class to compute rapially averaged structure factor of a 2D image""" 
    def __init__(self, shape): 
     assert len(shape) == 2, "only 2D images implemented" 
     L = max(shape) 
     self.qs = np.fft.fftfreq(L)[:L/2] 
     self.dists = np.sqrt(np.fft.fftfreq(shape[0])[:,None]**2 + np.fft.fftfreq(shape[1])**2) 
     self.dcount = np.histogram(self.dists.ravel(), bins=self.qs)[0] 
     self.has_window = False 

    def set_window(self,w='hanning'): 
     if w == False: 
      self.has_window = False 
     elif hasattr(np, w): 
      self.window = [getattr(np,w)(self.dists.shape[0])[:,None], getattr(np,w)(self.dists.shape[1])] 
      self.has_window = True 
     elif isinstance(w, np.ndarray): 
      assert w.shape == self.dists.shape 
      self.window = w 
      self.has_window = True 

    def windowing(self, im): 
     if not self.has_window: 
      return im 
     if isinstance(self.window, np.ndarray): 
      return numexpr.evaluate('im*w', {'im':im, 'w':self.window}) 
     return numexpr.evaluate('im*w0*w1', {'im':im, 'w0':self.window[0], 'w1':self.window[1]}) 

    def __call__(self, im): 
     spectrum = numexpr.evaluate(
      'real(abs(f))**2', 
      {'f':np.fft.fft2(self.windowing(im))} 
      ) 
     return np.histogram(self.dists.ravel(), bins=self.qs, weights=spectrum.ravel())[0]/self.dcount 

После этого импортируется, ваша проблема пишет:

size_patch=32 
isf = ImageStructureFactor((size_patch,size_patch)) 
psd1D_obs = np.zeroes((len(patches), len(isf.qs)-1)) 
for i, patch in enumerate(patches): 
    psd1D_obs[i] = isf(patch) 

где я предположил, что у вас есть список патчей. Если вы хотите, кадрирование, просто добавьте

isf.set_window() 

после построения isf. Вы можете указать имя окна (см. numpy.windowing) в виде строки.

Частоты доступны как isf.qs. Согласно руке numpy, первая не нулевая частота равна 1/size_patch, что, скорее, омега пульсации, чем частота f. Чтобы получить (пространственную) частоту, вам нужно умножить на 2 * pi.

Answer 2

np.fft.fftfreq(N, d = spacing) возвращает частоты выборок в циклах/интервалах. Если вы хотите, чтобы угловое волновое число просто умножалось на 2 * np.pi.

Вы также, скорее всего, захотите усреднить угловое значение при уменьшении 2d fft для 1-го представления. Если вы хотите использовать функцию удобства, которая красиво обертывает все эти детали, взгляните на https://github.com/keflavich/agpy/blob/master/AG_fft_tools/psds.py