# Simulation Box Power Spectrum/Multipoles (`FFTPower`

)¶

The `FFTPower`

class computes the 1d power spectrum \(P(k)\), 2d
power spectrum \(P(k,\mu)\), and/or multipoles \(P_\ell(k)\) for data
in a simulation box, using a Fast Fourier Transform (FFT). Here, we provide
a brief overview of the algorithm itself as well as the key things to know for
the user to get up and running quickly.

Note

To jump right into the `FFTPower`

algorithm, see this
cookbook recipe for a detailed
walk-through of the `FFTPower`

algorithm.

## The Algorithm¶

The steps involved in computing the power spectrum via `FFTPower`

are as follows:

**Generate data on a mesh**Data must be painted on to a discrete mesh to compute the power spectrum. There are several ways to generate data on a mesh (see Creating a Mesh), but the most common is painting a discrete catalog of objects on to a mesh (see Converting a CatalogSource to a Mesh and Painting Catalogs to a Mesh). The

`FFTPower`

class accepts input data in either the form of a`MeshSource`

or a`CatalogSource`

. In the latter case, the catalog is automatically converted to a mesh using the default parameters of the`to_mesh()`

function.When converting from a catalog to a mesh, users can customize the painting procedure via the options of the

`to_mesh()`

function. These options have important effects on the resulting power spectrum of the field in Fourier space. See Converting a CatalogSource to a Mesh for more details.**FFT the mesh to Fourier space**Once the density field is painted to the mesh, the Fourier transform of the field \(\delta(\vx)\) is performed in parallel to obtain the complex modes of the overdensity field, \(\delta(\vk)\). The field is stored using the

`ComplexField`

object.**Generate the 3D power spectrum on the mesh**The 3D power spectrum field is computed on the mesh, using

\[P(\mathbf{k}) = \delta(\mathbf{k}) \cdot \delta^\star(\mathbf{k}),\]where \(\delta^\star (\mathbf{k})\) is the complex conjugate of \(\delta(\mathbf{k})\).

**Perform the binning in the specified basis**Finally, the 3D power defined on the mesh \(P(\mathbf{k})\) is binned using the basis specified by the user. The available options for binning are:

1D binning as a function of wavenumber \(k\)

2D binning as a function of wavenumber \(k\) and cosine of the angle to the line-of-sight \(\mu\)

Multipole binning as a function of \(k\) and multipole number \(\ell\)

## The Functionality¶

Users can compute various quantities using the `FFTPower`

. We’ll discuss
the available functionality briefly in the sub-sections below.

### Auto power spectra and cross power spectra¶

Both auto and cross spectra are supported. Users can compute cross power spectra
by passing a second mesh object to the `FFTPower`

class using
the `second`

keyword. The first mesh object should always be specified as
the `first`

argument.

### 1D Power Spectrum, \(P(k)\)¶

The 1D power spectrum \(P(k)\) can be computed by specifying the
`mode`

argument as “1d”. The wavenumber binning will be linear, and can be
customized by specifying the `dk`

and `kmin`

attributes. By default,
the edge of the last wavenumber bin is the
Nyquist frequency, given
by \(k_\mathrm{Nyq} = \pi N_\mathrm{mesh} / L_\mathrm{box}\). If `dk`

is not specified, then the fundamental mode of the box is used:
\(2\pi/L_\mathrm{box}\).

### 2D Power Spectrum, \(P(k,\mu)\)¶

The 2D power spectrum \(P(k,\mu)\) can be computed by specifying the
`mode`

argument as “2d”. The number of \(\mu\) bins is specified via
the `Nmu`

keyword. The bins range from \(\mu=0\) to \(\mu=1\).

### Multipoles of \(P(k,\mu)\)¶

The `FFTPower`

class can also compute the multipoles of the 2D power
spectrum, defined as

where \(\mathcal{L}_\ell\) is the Legendre polynomial of order
\(\ell\). Users can specify which multipoles they wish to compute
by passing a list of the desired \(\ell\) values as the `poles`

keyword to the `FFTPower`

class.

For example, we can compute both \(P(k,\mu)\) and \(P_\ell(k)\) for a uniform catalog of objects using:

```
[2]:
```

```
from nbodykit.lab import UniformCatalog, FFTPower
cat = UniformCatalog(nbar=100, BoxSize=1.0, seed=42)
r = FFTPower(cat, mode='2d', Nmesh=32, Nmu=5, poles=[0,2,4])
```

```
/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.
from ._conv import register_converters as _register_converters
```

## The Results¶

The power spectrum results are stored in two attributes of the
initialized `FFTPower`

object:
`power`

and `poles`

. These attributes are
`BinnedStatistic`

objects, which
behave like structured numpy arrays and store
the measured results on a coordinate grid defined by the bins.
See Analyzing your Results for a full tutorial on using
the `BinnedStatistic`

class.

The `power`

attribute stores the following variables:

- k :
the mean value for each

`k`

bin

- muif
`mode=2d`

the mean value for each

`mu`

bin

- muif
- power :
complex array storing the real and imaginary components of the power

- modes :
the number of Fourier modes averaged together in each bin

The `poles`

attribute stores the following variables:

- k :
the mean value for each

`k`

bin

- power_L :
complex array storing the real and imaginary components for the \(\ell=L\) multipole

- modes :
the number of Fourier modes averaged together in each bin

Note that measured power results for bins where `modes`

is zero (no data points
to average over) are set to `NaN`

.

In our example, the `power`

and `poles`

attributes are:

```
[3]:
```

```
# the 2D power spectrum results
print("power = ", r.power)
print("variables = ", r.power.variables)
for name in r.power.variables:
var = r.power[name]
print("'%s' has shape %s and dtype %s" %(name, var.shape, var.dtype))
```

```
power = <BinnedStatistic: dims: (k: 16, mu: 5), variables: ('k', 'mu', 'power', 'modes')>
variables = ['k', 'mu', 'power', 'modes']
'k' has shape (16, 5) and dtype float64
'mu' has shape (16, 5) and dtype float64
'power' has shape (16, 5) and dtype complex128
'modes' has shape (16, 5) and dtype int64
```

```
[4]:
```

```
# the multipole results
print("poles = ", r.poles)
print("variables = ", r.poles.variables)
for name in r.poles.variables:
var = r.poles[name]
print("'%s' has shape %s and dtype %s" %(name, var.shape, var.dtype))
```

```
poles = <BinnedStatistic: dims: (k: 16), variables: 5 total>
variables = ['k', 'power_0', 'power_2', 'power_4', 'modes']
'k' has shape (16,) and dtype float64
'power_0' has shape (16,) and dtype complex128
'power_2' has shape (16,) and dtype complex128
'power_4' has shape (16,) and dtype complex128
'modes' has shape (16,) and dtype int64
```

These attributes also store meta-data computed during the power calculation
in the `attrs`

dictionary. Most importantly, the `shotnoise`

key
gives the Poisson shot noise, \(P_\mathrm{shot} = V / N\), where *V*
is the volume of the simulation box and *N* is the number of objects. The keys
`N1`

and `N2`

store the number of objects.

In our example, the meta-data is:

```
[5]:
```

```
for k in r.power.attrs:
print("%s = %s" %(k, str(r.power.attrs[k])))
```

```
Nmesh = [32 32 32]
BoxSize = [1. 1. 1.]
Lx = 1.0
Ly = 1.0
Lz = 1.0
volume = 1.0
mode = 2d
los = [0, 0, 1]
Nmu = 5
poles = [0, 2, 4]
dk = 6.283185307179586
kmin = 0.0
N1 = 96
N2 = 96
shotnoise = 0.010416666666666666
```

Note

The shot noise is not subtracted from any measured results. Users can
access the Poisson shot noise value in the meta-data `attrs`

dictionary.

## Saving and Loading¶

Results can easily be saved and loaded from disk in a reproducible manner
using the `FFTPower.save()`

and `FFTPower.load()`

functions.
The `save`

function stores the state of the algorithm,
including the meta-data in the `FFTPower.attrs`

dictionary, in a
JSON plaintext format.

```
[6]:
```

```
# save to file
r.save("fftpower-example.json")
# load from file
r2 = FFTPower.load("fftpower-example.json")
print("power = ", r2.power)
print("poles = ", r2.poles)
print("attrs = ", r2.attrs)
```

```
power = <BinnedStatistic: dims: (k: 16, mu: 5), variables: ('k', 'mu', 'power', 'modes')>
poles = <BinnedStatistic: dims: (k: 16), variables: 5 total>
attrs = {'Nmesh': array([32, 32, 32]), 'BoxSize': array([1., 1., 1.]), 'Lx': 1.0, 'Ly': 1.0, 'Lz': 1.0, 'volume': 1.0, 'mode': '2d', 'los': [0, 0, 1], 'Nmu': 5, 'poles': [0, 2, 4], 'dk': 6.283185307179586, 'kmin': 0.0, 'N1': 96, 'N2': 96, 'shotnoise': 0.010416666666666666}
```

## Common Pitfalls¶

The default configuration of nbodykit should lead to reasonable results
when using the `FFTPower`

algorithm. When performing custom, more complex
analyses, some of the more common pitfalls are:

When the results of

`FFTPower`

do not seem to make sense, the most common culprit is usually the configuration of the mesh, and whether or not the mesh is “compensated”. In the language of nbodykit, “compensated” refers to whether the effects of the interpolation window used to paint the density field have been de-convolved in Fourier space. See the Converting a CatalogSource to a Mesh section for detailed notes on this procedure.Be wary of normalization issues when painting weighted density fields. See Painting Catalogs to a Mesh for further details on painting meshes and Applying Functions to the Mesh for notes on applying arbitrary functions to the mesh while painting. See this cookbook recipe for examples of more advanced painting uses.