I want to plot the power spectrum of natural images so I've been using code from this resource: https://bertvandenbroucke.netlify.app/2019/05/24/computing-a-power-spectrum-in-python/ but I just want to verify that this is the correct implementation. Here's my specific code:
def power_spectrum(im_path,plot=True,transform=True):
if type(im_path) == str:
image = PIL.Image.open(im_path)
if transform:
image = torchvision.transforms.functional.center_crop(image,(256,256))
image = np.array(torchvision.transforms.functional.rgb_to_grayscale(image))
# print(image.shape)
else:
image = np.array(image)
else:
image = im_path
# print(np.array(image).shape)
npix = image.shape[0]
fourier_image = np.fft.fftn(image)
fourier_amplitudes = np.abs(fourier_image)**2
kfreq = np.fft.fftfreq(npix) * npix
kfreq2D = np.meshgrid(kfreq, kfreq)
knrm = np.sqrt(kfreq2D[0]**2 + kfreq2D[1]**2)
knrm = knrm.flatten()
fourier_amplitudes = fourier_amplitudes.flatten()
kbins = np.arange(0.5, npix//2+1, 1.)
kvals = 0.5 * (kbins[1:] + kbins[:-1])
Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes,
statistic = "mean",
bins = kbins)
Abins *= np.pi * (kbins[1:]**2 - kbins[:-1]**2)
if plot:
plot_power_spectra([image],[kvals],[Abins])
return kvals,Abins
And here is an example result:
