mascado.distortions.powerspectrum module

mascado.distortions.powerspectrum.psd_histogram(spectrum)

Bin all pixels of spectrum with same frequency.

The pixels on the border of a square centered on the spectrum’s center are considered to have the same frequency. Then Nyquist frequency corresponds to all outermost pixels.

These pixels of the same frequency are summed up into one bin.

Parameters

spectrum ((N, N)-shaped float array) – Power spectrum

Returns

bins – Binned spectrum.

Return type

(N//2+1)-shaped float array

mascado.distortions.powerspectrum.psd_histogram_cumulative(spectrum)

Bin all pixels of spectrum up to some frequency.

Cumulative version of psd_histogram. The pixels inside and on the border of a square are summed up into one bin.

Parameters

spectrum ((N, N)-shaped float array) – Power spectrum

Returns

bins – Binned spectrum.

Return type

(N//2+1)-shaped float array

mascado.distortions.powerspectrum.psd_spectrums(vf, posscale, grid=None)

Calculate power spectrum in both directions of vector field.

Parameters

• vf (mascado.distortions.polynomials.PolyVectorField) – 2D vector field.

• posscale (float) – Scale of vector field from [-1, 1] to arcsec. 2*posscale = FOV.

• grid (((N, N), (N, N)) 2-tuple of square float arrays) – Regular quadratic meshgrid (xx, yy) with ij-indexing covering \([-1,1]^2\) to evaluate vector field on. If no grid is given, one is created with side length of vector field degree times 4 + 1.

Returns

• xpsd ((N, N) float array)

• ypsd ((N, N) float array) – Power spectra of x- and y-component of vector field. The zeroth frequency / offset is in the center. Every pixel to the left, right, top, or bottom are increasing in frequency by \(1/\text{FOV}\).

Notes

The power spectrum \(P\) in spatial frequencies \(k\) of the distortion pattern \(D\) is given by

\[P[D](k) = \left| \mathcal F[D](k) \right|^2\]

Where the Fourier transform is

\[\mathcal F[D](k) = \int_{-\text{FOV}/2}^{\text{FOV}/2} D(x)\, e^{2 \pi i k x} \,\mathrm d x\]

With the substitution \(x^\prime=\frac{x}{\text{FOV}/2}\) we get \(D\) in normalized coordinates

\[\mathcal F[D](k) = \frac{1}{2}\text{FOV} \int_{-1}^1 \underbrace{D(x^\prime\,\text{FOV}/2)}_{D_\text{norm}(x^\prime)}\, e^{2 \pi i k x^\prime\, \text{FOV}/2} \,\mathrm d x^\prime\]

The input to NumPy’s discrete Fourier transform is an image and it corresponds to integration bounds from 0 to 1, so another substitution \(x^{\prime\prime}=\frac{x^\prime+1}{2}\) is needed:

\[\begin{split}\mathcal F[D](k) &= \text{FOV} \int_0^1 \underbrace{D(\text{FOV}\, x^{\prime\prime} - \text{FOV}/2)} _{D_\text{norm}^\text{img}(x^{\prime\prime})}\, e^{2 \pi i k (\text{FOV}\, x^{\prime\prime} - \text{FOV}/2)} \,\mathrm d x^{\prime\prime} \\[1em] &= \text{FOV} \cdot \mathcal F [D_\text{norm}^\text{img}](\text{FOV}\cdot k) \cdot \underbrace{e^{-2 \pi i k\, \text{FOV}/2}} _{\text{neglect phase}}\end{split}\]

Here we see, that the \(n\)-th frequency in the transformed result is the spatial frequency in units of \(n/\text{FOV}\).