Plot power

[MMesch]

This is a proposition to provide a single function to plot power spectra.

• new function .plot_powerspectrum(unit=xxx) (i.e. summed power)
  • .plot_powerspectrum(unit='per_l') plots the standard total power of a degree l
  • .plot_powerspectrum(unit='per_lm') plots the power of the individual coefficients in a 2d plot (see image below)
  • .plot_powerspectrum(unit='per_band', bandwidth=2) plots the summed power in an octave band

I propose to keep all sums it in coefficient normalization as opposed to always converting to 4pi.

The next step is to define powerspectraldensity in a clean way, although I think this will be quite difficult because 1) density implies always a 'power per something' and 2) it implies some kind of averaging (in frequency or in space) over something. E.g. right now SHPowerSpectrumDensity is the power per lm averaged over one particular l. Many other power spectral densities are possible. In the random model world, the global power per degree l (and a small neighbourhood) can be seen as the power per degree l averaged over local subregions of the model. I don't know what the best clean distinction is.

new plots:

Earth spectrum (linear plot). Slight North-South anisotropy at degrees 3-5

This is the Mars spectrum (log plot), zonal degree 1/2...

[MMesch]

before merging, I would like to implement the grid search for symmetry axes as well.

[MMesch]

anisotropy spectrum (Mars example)

[MMesch]

And the final addition: grid search for symmetry axes (blue = zonal, red = tesseral)

[MMesch]

https://github.com/MMesch/SHTOOLS/blob/plot_power/examples/notebooks/vesta_topography.ipynb

[MarkWieczorek]

• Can we combine all the plot_powerspectrum(), plot_power_perdegree(), plot_powerperlm(), and plot_powerperband() into one function, similar to get_powerspectrum()?
• Should get_powerspectrum() also have an option for per_band? Then we can delete the function get_powerperband().
• This description for the power spectrum is not entirely accurate (line 822): the sum_m abs(c_lm)**2 at degree l is computed. This is the power per degree, but only when using 4pi normalization. _get_powerperdegree() correctly calculates the total power for each degree by taking into account the normalization used for the coefficients (i.e., if you look at the function, you will see that schmidt and orth divide by (2l+1) and 4pi.)
• line 1047: def return_coeffs(self, normalization=None, csphase=None). normalization and csphase have to be set for the classes to work properly. Either we force the user to input these two everytime, or we provide default values, which were 4pi and nocsphase. The documentation said in multiple places that the default is 4pi and no phase, so my preference would be to set these two values to the default values. The problem is that the average user is not going to know what to use for csphase.
• For get_anisotropyspectrum() and get_symmetries(), and plot_symmetries() there is not enough documentation for me to understand what these do, and what they are calculating. Can you provide some references? My immediate reaction was that these might be a bit too esoteric, but if people use these in the literature, we should provide them.

[MMesch]

• Can we combine all the plot_powerspectrum(), plot_power_perdegree(), plot_powerperlm(), and plot_powerperband() into one function, similar to get_powerspectrum()?

yes. That should be the goal. Ultimately the underlying functions for particular spectra could become private.

• Should get_powerspectrum() also have an option for per_band? Then we can delete the function
  get_powerperband().

yes

• This description for the power spectrum is not entirely accurate (line 822): the sum_m abs(c_lm)**2 at degree l is computed. This is the power per degree, but only when using 4pi normalization. _get_powerperdegree() correctly calculates the total power for each degree by taking into account the normalization used for the coefficients (i.e., if you look at the function, you will see that schmidt and orth divide by (2l+1) and 4pi.)

this is a question of definition. In my opinion the power spectrum should (!) depend on the normalization. It merely states how much energy you have in different basis functions. If the basis functions are normalized differently, the power spectrum needs to change as well. I agree that the 4pi normalized power is probably easiest to interpret physically because it corresponds to the harmonics variance in space. The unit normalized power corresponds to the total energy of the harmonic in space. The other normalizations I don't know about.

• line 1047: def return_coeffs(self, normalization=None, csphase=None). normalization and csphase have to be set for the classes to work properly. Either we force the user to input these two everytime, or we provide default values, which were 4pi and nocsphase. The documentation said in multiple places that the default is 4pi and no phase, so my preference would be to set these two values to the default values. The problem is that the average user is not going to know what to use for csphase.

In my opinion, return_coeffs should return the same normalization as the coefficients have. This means that a call to return_coeffs without arguments initializes a clean copy of SHCoeffs (much better than using the copy module).

• For get_anisotropyspectrum() and get_symmetries(), and plot_symmetries() there is not enough documentation for me to understand what these do, and what they are calculating. Can you provide some references? My immediate reaction was that these might be a bit too esoteric, but if people use these in the literature, we should provide them.

OK. I'll provide better documentation. These methods are used e.g. in cosmology. Everything they do is to analyse whether energy is concentrated in high m (sectorial harmonics) or low m (zonal) coefficients. I think it is a nice addition that allows to go past the usual per degree analysis that throws away all information about orientation that is contained in a 2D spectrum.

[MarkWieczorek]

For the definition of the power spectrum, what I have been using is simply Parsival's theorem:

1/(4 pi) Integral f**2 dOmega = sum S(l)

It is S(l) that I call the power spectrum. With this definition, the numerical value of S(l) is independent of normalization. However, the exact equation that relates S to the coefficients does depend on the normalization.

We could define an (internal?) variable that is the L2 norm of the vector (i.e., sum of the coefficients squared) and then relate this to S(l).

[MMesch]

I think that's a good solution. Let's have a flag that states switches on the 4pi normalization or not instead of different subroutines. The 4pi normalized power has units power per unit area and per l, whereas the unit power has units (total) power per l.

Parseval's theorem can simply be stated differently without the 1/(4 pi) and you will have a different definition ... Many (!) papers use the unit normalized L2 norm of the coefficients as 'power spectrum' (e.g. as in Dahlen & Tromp 1998). Therefore I don't think that the 4pi is a general rule .

On the web I found this:

Wikipedia has a very unclear definition in my opinion:

When the energy of the signal is concentrated around a finite time interval, especially if its total energy is finite, one may compute the energy spectral density. More commonly used is the power spectral density (or simply power spectrum)

From Wolfram Mathworld it sounds like the power spectrum should be per unit area, i.e. 4pi normalized. Anyway that's just a choice of units:

For a given signal, the power spectrum gives a plot of the portion of a signal's power (energy per unit time) falling within given frequency bins. The most common way of generating a power spectrum is by using a discrete Fourier transform, but other techniques such as the maximum entropy method can also be used.

[MarkWieczorek]

Here is where the 4pi comes from

Power = [ Integral f**2  r**2 sin(theta) dthete dphi ] / [4pi r**2] = sum S(l)

A lot of people work on the "unit sphere" where r=1, and in this case it is easy to find an argument to ignore the 4pi factor. However, as soon as the sphere has a radius that is not 1, it becomes more difficult to ignore this.

[MMesch]

You normalize power by the area of the sphere in this equation. This is an arbitrary definition that many people chose differently. They will have different units for their power spectrum. In my opinion this is just a matter of units and defintiions. We cannot claim that 4pi is what has to be used, even though it has really sensible units.

EDIT: quick example. Think about quantum mechanics, where probability densities are expressed in spherical harmonics. In this case it makes really sense to use unit normalized spherical harmonics and the L2 norm of their coefficients as power, because it is the total power on the sphere that is interesting and not necessarily the power per unit area.

[MarkWieczorek]

Can you write the equation for the power that you are proposing for the function f(theta, phi), without any reference to spherical harmonics? Here are my two definitions for power and energy:

Power = Integral f(theta, phi)**2 r**2 sin(theta) dtheta dpi / (4pi r**2)
Energy = Integral f(theta, phi)**2 r**2 sin(theta) dtheta dpi

[MMesch]

Well I have seen both definitions called 'power'. This is my intuition (now):

Energy = Integral f(theta, phi)**2 r**2 sin(theta) dtheta dpi
Energy per unit area (Power=Variance) = Integral f(theta, phi)**2 r**2 sin(theta) dtheta dpi / (4pi r**2)
Energy per unit area and degree l = Integral f_l(theta, phi)**2 r**2 sin(theta) dtheta dpi / (4pi r**2)
Energy per unit area and degree l and order m = Integral f_lm(theta, phi)**2 r**2 sin(theta) dtheta dpi / (4pi r**2)

I think you are right Mark, at least following the standard definition from physics:

power is the rate of doing work. It is the amount of energy consumed per unit time [i.e. area in our case]

Maybe we should just call the proposed function energy then?

Dahlen & Tromp 1998 has yet another definition or is just incorrect (?) Eq. B95 and following.
The L2 norm of the orthonormal coefficients is defined as:

variance or power per degree l and per unit area of the function Psi

This doesn't make sense in my opinion.

[MMesch]

another proposal:

.get_spectrum('energy_per_l')
.get_spectrum('energy_per_logl')
.get_spectrum('energy_per_lm')
.get_spectrum('power_per_l')
.get_spectrum('power_per_lm')
.get_spectrum('power_per_logl')

what do you think?

[MMesch]

here one more definition of power, as coefficient variance of orthonormal harmonics:
http://ada7.cosmostat.org/ADA7_Tutorial_6.pdf
http://healpix.jpl.nasa.gov/pdf/intro.pdf

often it is just called spectral density function which avoids the physical connotations of power and energy... E.g.:
http://publications.nr.no/directdownload/rask/old/917_Rapport.pdf

To sum up, there are different definitions in different communities. SHTOOLS seems to follow the signal processing definitions.

[MarkWieczorek]

For the publication that we will eventually write, we will have to summarize all of these definitions. The best way to do this will probably be to provide the equation that relates the quantity measured in the space domain to the corresponding equations using spherical harmonic coefficents. What makes this so confusing is that most people only provide the equation for
S(l) ~ sum(coefficents**2)
but not the relation between
sum S(l) ~ integral (function**2)

I wonder if we should break this pull request into several smaller ones:

• return_coeffs : I think that I am ok with this, but would like to verify that this is treated consistently in all routines.
• get_powerspectrum, plot_powerspectrum : This is the heart of the matter.
• get_anisotropyspectrum() and get_symmetries(), and plot_symmetries(). This will probably be the last thing to be merged, and if we split this, that should speed things up.
• a page somewhere (in the PR comments?) summarizing all the different definitions so that we have a record of this later.

[MMesch]

I'll split it up. Pls leave this PR open ...

[MarkWieczorek]

Can you refactor this pull request, given that most of the powerspectrum stuff has been added to master?

[MarkWieczorek]

I'm going to close this issue as no progress has been made in over a year, and many aspects of this have already been implemented. Concerning the anisotropy spectrum, I think that we might want to also look at the approach in Bills and Lemoine 1995 who define this as the ratio of the RMS slopes in the NS and EW direction.