CODEWITHSUNDEEP
numpyfftfrequencydft
user3722235
user3722235

Reputation: 567

Frequency axis in a Numpy fft

I have a 32x32 pixels image that I would like to compute the 1D power spectrum (averaged along the y axis of the image). Here's what I do.

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

I have some trouble grasping what is precisely plotted as an x-axis in the power spectrum. Is it the wavenumber or the spatial frequency? What is the convention adopted here? Is the natural unit cycles/pixel? The doc for numpy.fft.fftfreq is a bit too vague to me. This might be a very simple question but I didn't find any clear answer anywhere. Could you please enlighten me?

Upvotes: 2

Views: 3905

Answers (2)

Andi
Andi

Reputation: 1273

np.fft.fftfreq(N, d = spacing) returns the sample frequencies in cycles / spacing. If you want to have the angular wavenumber instead simply multiply by 2 * np.pi.

You also most likely want to average angularly when reducing your 2d fft for a 1d representation. If you want to use a convenience function that nicely wraps all of these details, then take a look at https://github.com/keflavich/agpy/blob/master/AG_fft_tools/psds.py

Upvotes: 2

Mathieu Leocmach
Mathieu Leocmach

Reputation: 61

If you want to average a 2D spectrum, you should do a radial average, not along an axis. A radial average between r and r+dr is the sum of the values of the pixels whoes distance r to the center is between r and dr, divided by the number of those pixels. This is feasible with a loop, but I'd rather use numpy.histogram twice.

Here is a class that does the job, with possibility of windowing. It is optimised for the common case where you have many images (or patches) of the same size to treat. I also use numexpr to speed up what can be.

In the following I name my spatial frequency 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

Once this is imported, your problem writes:

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)

where I supposed that you had a list of patches. If you want windowing, just add

isf.set_window()

after constructing isf. You can specify the name of the window (see numpy.windowing) as a string.

The frequencies are available as isf.qs. According to numpy manual, the first non null frequency is 1/size_patch, thus rather a pulsation omega than a frequency f. To get a (spatial) frequency you need to multiply by 2*pi.

Upvotes: 0

Related Questions