Numpy fft中的频率轴答案

【问题标题】：Frequency axis in a Numpy fftNumpy fft中的频率轴
【发布时间】：2015-04-09 19:59:56
【问题描述】：

我有一个 32x32 像素的图像，我想计算一维功率谱（沿图像的 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

我在掌握功率谱中精确绘制为 x 轴的内容时遇到了一些麻烦。是波数还是空间频率？这里采用的约定是什么？是自然单位周期/像素吗？ numpy.fft.fftfreq 的文档对我来说有点模糊。这可能是一个非常简单的问题，但我在任何地方都没有找到任何明确的答案。可以请教一下吗？

【问题讨论】：

来自您的链接：“返回的浮点数组 f 包含以每单位样本间隔为周期的频率 bin 中心（开头为零）。”在您的情况下，这转化为每像素周期的空间频率。

标签： numpy fft frequency dft

【解决方案1】：

np.fft.fftfreq(N, d = spacing) 以周期/间隔返回采样频率。如果你想得到角波数，只需乘以2 * np.pi。

在为 1d 表示减少 2d fft 时，您也很可能希望平均角度。如果你想使用一个方便的函数来很好地包装所有这些细节，那么看看https://github.com/keflavich/agpy/blob/master/AG_fft_tools/psds.py

【讨论】：

【解决方案2】：

如果您想平均二维光谱，您应该进行径向平均，而不是沿轴。 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。

【讨论】：