## Computing a power spectrum in Python

A power spectrum is an analysis tool that is very often used to do a statistical analysis for a large, seemingly chaotic data set. Within astronomy, they are used in practically every field: the power spectrum of the cosmic microwave background radiation contains vital clues about the Big Bang and the history of the Universe, the power spectrum of the radio emission from galaxies can be used to characterise the distribution of interstellar gas within those galaxies, and the power spectrum of measured gas velocities in star-forming molecular clouds can be used to investigate turbulence within those clouds.

Yet, despite their widespread use, power spectra are not as easy to use as people think. First of all, there seems to be some confusion about what a power spectrum is, which I want to clear out below. Second of all, it is not as easy to compute the power spectrum for a given data set. IDL (a proprietary computing tool that is nonetheless often used in astronomy) has library functions that do it for you, but Python (a non-proprietary computing tool that is much wider used than IDL and has my personal preference) does not. In this post, I will explain how to compute a power spectrum using Python.

# What is a power spectrum?

But letâ€™s start with what a power spectrum actually is. Imagine you have an image of a set of clouds, like the one below (in grey scale for reasons that will become obvious later). The image represents a spatial map of cloud intensity: where the image is white, the cloud is strongly present, while it is less present in darker regions.

Clearly, these clouds have a very complicated, seemingly chaotic structure. They are made up of many small parts and have a fractal-like appearance: when you zoom in on a portion of the cloud near the edge, the zoomed in version of the image would look very similar to the original image.

Despite the chaos, we could still try to get some statistical information about the clouds in this image. We could for example try to determine what the typical size of a smaller cloud within the larger cloud is. To do this, we could e.g. decompose the image into different size scales. We can use some image manipulation tool to rescale the image to a very coarse version, like a 2x2 pixel image. This image will have pixels that contain the dominant colour within each large pixel in the original image. If we subtract this dominant colour from the original image, we are left with a version of the image without these largest structures. We can then repeat this procedure for increasingly more refined versions of the original image (always subtracting out the dominant components in the coarser levels), until we end up at the original image resolution. Every image in between will now contain a measurement of how strong features in the original image are at that specific resolution, and the size scale associated with that resolution tells us what the size scale of those features is.

A power spectrum is in essence a measure of the strength of the different features at the different resolutions. For each resolution, it is defined as the *variance* of the features on that resolution, i.e. how much total spread there is in values. A large spread means a strong signal and hence a lot of variation on this scale.

In order to do the proposed decomposition in resolutions, we require a mathematical technique called a Fourier transform. The signal on each resolution is measured by setting up a periodic function consisting of sines and cosines, with the period of these sines and cosines, called , linked to the resolution (a period corresponds to a sine wave that fits exactly once inside the frame of the image). For each value of , we get two *Fourier amplitudes*: the amplitudes of the corresponding sine and cosine waves that best fit the pixels in one image dimension. Since the image has two dimensions, itself has two components, so that we end up with four Fourier amplitudes per value of . The sine and cosine amplitudes are usually combined into a single *complex* number for convenience.

In principle, this Fourier decomposition can be done for an infinite number of values. In practice, there is however a clear limit: when the period of the sine waves is equal to the number of pixels in one dimension, then each oscillation of the sine wave exactly covers one pixel. If we where to increase the resolution even more, then multiple sine waves would cover a single pixel, and the resulting amplitude would no longer change by increasing the resolution even more. The maximum frequency that is useful in an image is hence , with the number of pixels in the image. In practice, this means we need exactly the same number of complex pixels in Fourier space to represent our image, as there are pixels in the original image.

For a multi-dimensional Fourier decomposition, there are many combinations of that represent the same resolution . However, they all encode the same information, so that we usually bin these together. The appropriate measure for the power for a specific value is than the average power within the corresponding bin, multiplied with the volume of that bin in space. For a two dimensional Fourier decomposition, this bin volume increases linearly with , in three dimensions it even increases as .

Note that a Fourier decomposition is not the only possible way of decomposing a data set. Other spatial decompositions are possible, and the power spectrum for those can also be called *the* power spectrum. When the data set is taken on the surface of the sphere (e.g. the cosmic microwave background data that is measured on the celestial sphere), then a decomposition in terms of *spherical harmonics* is used. Again, it is possible to calculate a power in spherical harmonics space that links to a specific *angular scale*, and the corresponding curve of power versus angular scale is also called *the* power spectrum. Hence the confusion I mentioned before.

# Calculating the power spectrum in Python

In order to calculate the power spectrum for a data set, we have to do the following:

- Convert the data set into a suitable data array with the correct spatial layout.
- Take the Fourier transform of the data array in the number of dimensions of the data array.
- Construct a corresponding array of wave vectors with the same layout as the Fourier amplitudes.
- Bin the amplitudes of the Fourier signal into bins and compute the total variance within each bin.

Below, I will illustrate these steps for the cloud image introduced above.

## Constructing a data array

This first step depends a lot on the data set under study. For the cloud image, we want to read in the pixel data for the image. We can do this using `matplotlib`

:

```
import matplotlib.image as mpimg
image = mpimg.imread("clouds.png")
```

This will simply create a 1000x1000 `numpy.array`

of floating point values that represents the grey scale level of each pixel. You can visualise this list using `imshow`

:

```
import pylab as pl
pl.imshow(image)
pl.show()
```

Since this function will map the pixel values to the default colour scale, the result will look quite interesting.

## Taking the Fourier transform

To take the Fourier transform of our two dimensional image data array, we will use `numpy.fft`

. This library has a number of functions to handle one, two and multi dimensional data arrays. We will use the multi dimensional function here, as that makes it possible to generalise the technique to three dimensional data sets quite easily:

```
import numpy as np
fourier_image = np.fft.fftn(image)
```

The Fourier image array now contains the *complex valued* amplitudes of all the Fourier components. We are only interested in the size of these amplitudes. We will further assume that the average amplitude is zero, so that we only require the square of the amplitudes to compute the variances. We can get these values like this

```
fourier_amplitudes = np.abs(fourier_image)**2
```

# Constructing a wave vector array

To bin the results found above in space, we need to know what the layout of the return value of `numpy.fft.fftn`

is: what is the wave vector corresponding to an element with indices `i`

and `j`

in the return array? If we were to make things complicated for ourselves, we could try to figure this out from the online documentation.

Alternatively, we could not worry about this, and just use the utility function provided by `numpy`

:

```
kfreq = np.fft.fftfreq(1000) * 1000
```

This will automatically return a one dimensional array containing the wave vectors for the `numpy.fft.fftn`

call, in the correct order. By default, the wave vectors are given as a fraction of 1, by multiplying with the total number of pixels, we convert them to a pixel frequency. To convert this to a two dimensional array matching the layout of the two dimensional Fourier image, we can use `numpy.meshgrid`

:

```
kfreq2D = np.meshgrid(kfreq, kfreq)
```

Finally, we are not really interested in the actual wave vectors, but rather in their norm:

```
knrm = np.sqrt(kfreq2D[0]**2 + kfreq2D[1]**2)
```

For what follows, we no longer need the wave vector norms or Fourier image to be laid out as a two dimensional array, so we will flatten them:

```
knrm = knrm.flatten()
fourier_amplitudes = fourier_amplitudes.flatten()
```

## Creating the power spectrum

To bin the amplitudes in space, we need to set up wave number bins. We will create integer value bins, as is common:

```
kbins = np.arange(0.5, 501., 1.)
```

Note that the maximum wave number will equal half the pixel size of the image. This is because half of the Fourier frequencies can be mapped back to negative wave numbers that have the same norm as their positive counterpart. The `kbin`

array will contain the start and end points of all bins; the corresponding values are the midpoints of these bins:

```
kvals = 0.5 * (kbins[1:] + kbins[:-1])
```

To compute the average Fourier amplitude (squared) in each bin, we can use `scipy.stats`

:

```
import scipy.stats as stats
Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes,
statistic = "mean",
bins = kbins)
```

Remember that we want the *total* variance within each bin. Right now, we only have the average power. To get the total power, we need to multiply with the volume in each bin:

```
Abins *= 4. * np.pi / 3. * (kbins[1:]**3 - kbins[:-1]**3)
```

Finally, we can plot the resulting power spectrum as a function of wave number, (typically plotted on a double logarithmic scale):

This power spectrum has some interesting features. First of all, we can see that most of the power is located at large scales (small wave numbers). This makes sense, as the image is dominated by the large cloud structures. Towards smaller scales (larger wave numbers), the power drops of almost linearly. This is an expression of the fractal nature of the cloud: at lower resolutions the same types of patterns return, but at increasingly lower signal. There is a huge spike in the power spectrum at a wave number of , or a size scale of about 3 pixels. This is the scale at which the smallest cloud features at the edges of the clouds manifest.

# Putting it all together

Below is the full script to plot the power spectrum for the cloud image. I hope the part that calculates the power spectrum might be useful for other applications. An extension to higher dimensions is straightforward; simply replace the `knrm`

calculation with a higher order equivalent.

```
import matplotlib.image as mpimg
import numpy as np
image = mpimg.imread("clouds.png")
fourier_image = np.fft.fftn(image)
fourier_amplitudes = np.abs(fourier_image)**2
kfreq = np.fft.fftfreq(1000) * 1000
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, 501., 1.)
kvals = 0.5 * (kbins[1:] + kbins[:-1])
Abins, _, _ = stats.binned_statistic(knrm, fourier_amplitudes,
statistic = "mean",
bins = kbins)
Abins *= 4. * np.pi / 3. * (kbins[1:]**3 - kbins[:-1]**3)
pl.loglog(kvals, Abins)
pl.xlabel("$k$")
pl.ylabel("$P(k)$")
pl.tight_layout()
pl.savefig("cloud_power_spectrum.png", dpi = 300, bbox_inches = "tight")
```

Professional astronomer.