2
$\begingroup$

Imagine satellite images, these are irregular sampled in X and Y direction 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 spectra (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 2D power spectrum is radially averaged and plotted. Here the 2D power scalesthe 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 xrange(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

fractal noise

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()

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

amp = num.empty_like(k_rad)
amp[r0] = k_rad[r0] ** -beta[0]
amp[r0] /= amp[r0].max()

amp[r1] = k_rad[r1] ** -beta[1]
amp[r1] /= amp[r1].max()/amp[r0].min()

amp[r2] = k_rad[r2] ** -beta[2]
amp[r2] /= amp[r2].max()/amp[r1].min()

amp[k_rad == 0.] = amp.max()

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

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

ampE_slice = shift(amp)[amp.shape[0]/2, :amp.shape[1]/2]
kE_slice = shift(k_rad)[amp.shape[0]/2, :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
pspec[k_rad == 0.] = 0.

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 xrange(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 xrange(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)))
$\endgroup$
  • 1
    $\begingroup$ 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? $\endgroup$ – Maximilian Matthé Jan 12 '17 at 19:25
  • $\begingroup$ Okay, let me try $\endgroup$ – n1nj4 Jan 13 '17 at 9:21
  • $\begingroup$ This looks like a really useful example - shame it is a little bit light on detail about whats going inside $\endgroup$ – 4Oh4 Apr 30 at 17:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.