Radial Average for Power Spectrum

[santisoler]

First try to implement the Radial Average of a Power Spectrum (or any real fft function) proposed in #224 .

It needs a Power Spectrum in a square grid with odd number of points per axis in order to ensure that the center of the zero frequency matches the center of the grid. The radial average is made by concentric rings of inner radii f-df/2 and outer radii f+df/2, where f is each frecuency and df is the fundamental frecuency. The first ring is a circle containing only the zero frecuency point. The next is a ring surrounding this circle, and so on.

Also, it calculates the standard deviation for each average:
rms = \sqrt{ \sum (x_i - \overline{x})*_2 / nx_ny }

Here is an example:

import numpy as np
import matplotlib.pyplot as plt
from fatiando.datasets import fetch_bouguer_alps_egm
from fatiando.gridder import load_surfer
from fatiando.gravmag import fourier

# Load Bouguer Anomaly Data
fetch_bouguer_alps_egm(fname="alps-bouguer.grd")
x, y, bg, shape = load_surfer("alps-bouguer.grd")

## Puts Data in 2d array grid
x = x.reshape(shape)
y = y.reshape(shape)
bg = bg.reshape(shape)

## Cut the data to a square grid
x, y, bg = x[0:181, 19:200], y[0:181, 19:200], bg[0:181, 19:200]
dx = abs(x[0][1] - x[0][0])

## Calculates the Power Spectrum
fx = np.fft.fftfreq(181, abs(x[0][1] - x[0][0]))
fx, fy = np.meshgrid(fx, fx)
fft_bg = np.fft.fft2(bg)
log_power_spectrum = np.log(abs(fft_bg)**2)

## Calculates the Radial Average
freqs, average, rms = fourier.radial_average(log_power_spectrum, dx)

## PLOTS
plt.axes(aspect='equal')
plt.contourf(x, y, bg, 150)
plt.colorbar(label='mGal')
plt.title("Bouguer Anomaly")
plt.xlabel("m")
plt.ylabel("m")
plt.show()

fx = np.fft.fftshift(fx)
fy = np.fft.fftshift(fy)
log_power_spectrum = np.fft.fftshift(log_power_spectrum)
plt.axes(aspect='equal')
plt.contourf(fx, fy, log_power_spectrum, 150)
plt.colorbar()
plt.title("Log Power Spectrum")
plt.show()

plt.errorbar(freqs, average, yerr=rms, fmt='.')
plt.show()

[leouieda]

@santis19 Thanks for this! I'll have a look at this later on.

For now, your master branch was out of sync with the master of the main repo fatiando/fatiando. That's why you have merge conflicts. Please update your master branch and then merge it into this one. I can help you with that if you don't know how.

The gravmag.fourier module was abandoned in favor of gravmag.transform. All functions were moved there and this one should go there as well. You can re-use the FFT auxiliary functions to get the frequencies. There are some tricky things regarding the coordinate axis orientation so it's better to not re-implement things.

[santisoler]

Yes @leouieda, I started using the git via terminal instead of the web GUI. I think I clone my own old dated fork. I think I have to fork it again and then clone it to my desktop, right?

[leouieda]

You don't have to fork again, you can update your fork instead. Try these
instructions https://help.github.com/articles/syncing-a-fork/ Let me know
if you have any problems.
On Aug 11, 2015 8:47 AM, "Santiago Soler" [email protected] wrote:

Yes @leouieda https://github.com/leouieda, I started using the git via
terminal instead of the web GUI. I think I clone my own old dated fork. I
think I have to fork it again and then clone it to my desktop, right?

—
Reply to this email directly or view it on GitHub
#227 (comment).

[santisoler]

@leouieda
Thanks, I finally updated my fork.

The thing about coordinate axis regards that fatiando uses x for Northing and y for Easting? Because the _fftfreqs function returns the meshgrid in the inverted order.

On the other hand, for me to use that function, I have to ask for the x and y arrays, which I think that are unnecessary for the radial average, even confusing and can lead to mistakes in case that you work just with frequencies and ffts. Perhaps I should modify the meshgrid line to inverse the order of the axis. But it would be only for keeping a default order between the codes, because for a radial average is the same that you put an xy coordinates or a yx system, due to reflection simmetry.