# Calculate 1D Power Spectrum from 2D Images

Imagine satellite images, these are irregularly sampled in the X and Y directions, and the shapes are, of course are oddly off. We now want to estimate a 1D power spectrum from the whole image to estimate the atmospheric noise.

What works is the 2D power spectrum, power2DMean. Through:

$$\bar P_{r} = \frac{4 \pi}{M \cdot N} \cdot \frac{1}{R}\sum_{r=0}^{R} P_{k_r, l_r}$$

Well now I have difficulties to scale the 1D spectrum.

To prove the signal processing techniques, we are working on synthetic data in the code:

1. Fractal noise is created on from a given power spectrum (fractal noise, left plot; power spectrum, dashed, center plot)
2. Same power spectrum is created for a 1D spectrum (dashed, right plot)
3. Fractal noise is given as a 2D power spectrum, that is radially averaged and plotted. Here the 2D power scales the model, but the scaling factor of the 1D integral is unclear to me.

Interesting sections in the code below are power1D and power2DMean.

def power1D(k, N=256):
""" Integrating power over the radii k
Where scaling does not work....

power_interp interpolated power spectrum surface, unscaled.
"""
theta = num.linspace(-num.pi, num.pi, N, False)
power = num.empty_like(k)
for i in range(k.size):
kE = num.sin(theta) * k[i]
kN = num.cos(theta) * k[i]
power[i] = num.median(power_interp.ev(kN, kE) * k[i])
return power


Whole piece for completeness:

import matplotlib.pyplot as plt
import numpy as num
from scipy import interpolate

shift = num.fft.fftshift
fig, ax = plt.subplots(1, 3)

'''
Create fractal noise
'''
nN, nE = 1024, 512
dE, dN = 101., 134.  # Arbitrary values for sampling in dx and dy
amplitude = 50.

rfield = num.random.rand(nN, nE)
spec = num.fft.fft2(rfield)

regime = num.array([.15, .60, 1.])
beta = num.array([5./3, 8./3, 2./3])
beta += 1.  # Betas are defined for 1D PowerSpec, increasing dimension

kE = num.fft.fftfreq(nE, dE)
kN = num.fft.fftfreq(nN, dN)

k = kN if kN.size < kE.size else kE
k = k[k > 0]
k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2)

k0 = 0
k1 = regime[0] * k.max()
k2 = regime[1] * k.max()

amp[r0] /= amp[r0].max()

amp[r1] /= amp[r1].max()/amp[r0].min()

amp[r2] /= amp[r2].max()/amp[r1].min()

amp *= amplitude**2
spec *= num.sqrt(amp)  # We come from powerspec!
noise = num.abs(num.fft.ifft2(spec))

ampN_slice = shift(amp)[:int(amp.shape[0]/2), int(amp.shape[1]/2)]

ampE_slice = shift(amp)[int(amp.shape[0]/2), :int(amp.shape[1]/2)]

ax[0].imshow(noise)
ax[0].set_title('Fractal Noise')

'''
Model spec for 1D (used for comparison)
'''
k = num.linspace(max(kE[kE > 0.].min(), kN[kN > 0.].min()),
max(kE.max(), kN.max()), 512)
k1d = k

r0 = num.logical_and(k1d >= k0, k1d < k1)
r1 = num.logical_and(k1d >= k1, k1d < k2)
r2 = k1d >= k2

beta1d = num.array([5./3, 8./3, 2./3])
amp1d = num.zeros_like(k1d)

amp1d[r0] = k[r0] ** -(beta1d[0])
amp1d[r0] /= amp1d[r0].max()

amp1d[r1] = k[r1] ** -(beta1d[1])
amp1d[r1] /= amp1d[r1].max()/amp1d[r0].min()

s2 = k ** -beta[1]
amp1d[r2] = k[r2] ** -(beta1d[2])
amp1d[r2] /= amp1d[r2].max()/amp1d[r1].min()

amp1d *= amplitude**2  # We are in the powerspec

'''
Noise analysis from random 2D spectrum
'''
spec = shift(num.fft.fft2(noise))
pspec = num.abs(spec)**2

kE = shift(num.fft.fftfreq(spec.shape[1], dE))
kN = shift(num.fft.fftfreq(spec.shape[0], dN))
k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2)

kE = shift(num.fft.fftfreq(spec.shape[1], dE))
kN = shift(num.fft.fftfreq(spec.shape[0], dN))

power_interp = interpolate.RectBivariateSpline(kN, kE, pspec)

def power2DMean(k, N=256):
""" Mean 2D Power works! """
theta = num.linspace(-num.pi, num.pi, N, False)
power = num.empty_like(k)
for i in range(k.size):
kE = num.sin(theta) * k[i]
kN = num.cos(theta) * k[i]
power[i] = num.median(power_interp.ev(kN, kE) * 4 * num.pi)
# Median is more stable than the mean here
return power / pspec.size

def power1D(k, N=256):
""" Here we need to normalize the power over the radius
But scaling does not work """
theta = num.linspace(-num.pi, num.pi, N, False)
power = num.empty_like(k)
for i in range(k.size):
kE = num.sin(theta) * k[i]
kN = num.cos(theta) * k[i]
power[i] = num.median(power_interp.ev(kN, kE) * k[i])
return power

ax[1].loglog(kE_slice, ampE_slice, ls='--', label='Model Spec')
ax[1].loglog(k, power2DMean(k), alpha=.6, label='Noise Spec')
ax[1].legend()
ax[1].set_title('2D Mean Spectrum')
# dist = 10**(num.log10(ampE_slice) - num.log10(power1D(kE_slice)))

ax[2].loglog(k1d, amp1d, ls='--', label='Model Spec')
ax[2].loglog(k1d, power1D(k1d), alpha=.6, label='Noise Spec')
ax[2].legend()
ax[2].set_title('1D Averaged Spectrum')
plt.show()

# dist = 10**(num.log10(amp1d) - num.log10(power1D(k1d)))


• Actually, it's a lot code, and it's not easy to follow. Can you explain more, what a 1D Power Spectrum of a 2D image is supposed to be? Commented Jan 12, 2017 at 19:25
• Okay, let me try Commented Jan 13, 2017 at 9:21
• This looks like a really useful example - shame it is a little bit light on detail about whats going inside
– 4Oh4
Commented Apr 30, 2019 at 17:11
• power[i] = num.median(power_interp.ev(kN, kE) * k[i] * 256/(4 * num.pi)) Commented Mar 29, 2023 at 18:43

I am just reading Quantitative characterization of surface topography using spectral analysis and recalled seeing your post today. The paper makes distinction between $$C^{iso}$$ (radially averaged) and $$C^{pseudo-1D}$$ (radially averaged but scaled by the factor of $$\frac{q}{\pi}$$, for reasons described in there), that looks like the issue you have. This is fig. 2(f) from the paper: