Tools for generating and manipulating power spectra

The mxlib PSD utilities are based on the variance normalization for PSDs. For a 1-sided 1-dimensional PSD that is

\[ \sigma^2 = \int_0^{f_{max}} PSD(f) df \]

or for a 2-sided 1-D PSD

\[ \sigma^2 = \int_{f_{min}}^{f_{max}} PSD(f) df. \]

More generally

\[ \sigma^2 = \int_{\vec{k}} PSD(\vec{k}) d\vec{k}. \]

Classes
class	mx::sigproc::averagePeriodogram< realT >
	Calculate the average periodogram of a time-series. More...

struct	mx::sigproc::psdVarMean< paramsT >
	Calculate the variance of the mean for a process given its PSD. More...

Modules
	Autoregressive Power Spectra

	PSD Filter
	Filtering with a PSD to generate correlated noise.

Function Documentation

◆ augment1SidedPSD()

template<typename vectorTout , typename vectorTin >

void mx::sigproc::augment1SidedPSD	(	vectorTout &	psdTwoSided,
		vectorTin &	psdOneSided,
		bool	addZeroFreq = `false`,
		typename vectorTin::value_type	scale = `0.5`
	)

Augment a 1-sided PSD to standard 2-sided FFT form.

Allocates psdTwoSided to hold a flipped copy of psdOneSided. Default assumes that psdOneSided[0] corresponds to 0 frequency, but this can be changed by setting zeroFreq to a non-zero value. In this case psdTwoSided[0] is set to 0, and the augmented psd is shifted by 1.

To illustrate, the bins are re-ordered as:

* {1,2,3,4,5} --> {0,1,2,3,4,5,-4,-3,-2,-1}
*

The output is scaled so that the total power remains the same. The 0-freq and Nyquist freq are not scaled.

Entries in psdOneSided are cast to the value_type of psdTwoSided, for instance to allow for conversion to complex type.

Test:: Verify scaling and normalization of augment1SidedPSD [test doc]

Parameters

[out]	psdTwoSided	on return contains the FFT storage order copy of psdOneSided.
[in]	psdOneSided	the one-sided PSD to augment
[in]	addZeroFreq	[optional] set to true if psdOneSided does not contain a zero frequency component.
[in]	scale	[optional] value to scale the input by when copying to the output. The default 0.5 re-normalizes for a 2-sided PSD.

Definition at line 817 of file psdUtils.hpp.

Referenced by mx::AO::analysis::clAOLinearPredictor< _realT >::calcCoefficients(), mx::AO::analysis::fourierTemporalPSD< _realT, aosysT >::intensityPSD(), SCENARIO(), SCENARIO(), and mx::AO::analysis::speckleAmpPSD().

◆ augment1SidedPSDFreq()

template<typename T >

void mx::sigproc::augment1SidedPSDFreq	(	std::vector< T > &	freqTwoSided,
		std::vector< T > &	freqOneSided
	)

Augment a 1-sided frequency scale to standard FFT form.

Allocates freqTwoSided to hold a flipped copy of freqOneSided. If freqOneSided[0] is not 0, freqTwoSided[0] is set to 0, and the augmented frequency scale is shifted by 1.

Example:

{1,2,3,4,5} --> {0,1,2,3,4,5,-4,-3,-2,-1}

Test:: Verify scaling and normalization of augment1SidedPSD [test doc]

Parameters

[out]	freqTwoSided	on return contains the FFT storage order copy of freqOneSided.
[in]	freqOneSided	the one-sided frequency scale to augment

Definition at line 876 of file psdUtils.hpp.

Referenced by mx::AO::analysis::fourierTemporalPSD< _realT, aosysT >::intensityPSD(), SCENARIO(), and SCENARIO().

◆ freq_sampling()

template<class realT >

realT mx::sigproc::freq_sampling	(	size_t	dim,
		realT	f_max
	)

Calculates the frequency sampling for a grid given maximum dimension and maximum frequency.

The freq_sampling is \( \Delta f = f_{max}/ (0.5*dim) \) where \( f_{max} = 1/(2\Delta t) \) is the maximum frequency and \( dim \) is the size of the grid.

Parameters

[in]	dim	is the size of the grid
[in]	f_max	is the maximum frequency of the grid

Returns: the sampling interval \( \delta f \)

Template Parameters

realT is the real floating point type used for calculations.

Definition at line 196 of file psdUtils.hpp.

Referenced by mx::sigproc::frequency_grid1D(), and mx::sigproc::frequencyGrid().

◆ frequency_grid1D()

template<typename eigenArr >

void mx::sigproc::frequency_grid1D	(	eigenArr &	vec,
		typename eigenArr::Scalar	dt,
		bool	inverse = `false`
	)

Create a 1-D frequency grid.

Parameters

[out]	vec	the pre-allocated Eigen-type 1xN or Nx1 array, on return contains the frequency grid
[in]	dt	the temporal sampling of the time series
[in]	inverse	[optional] if true

Template Parameters

eigenArr the Eigen-like array type

Definition at line 211 of file psdUtils.hpp.

References mx::sigproc::freq_sampling().

◆ frequencyGrid() [1/4]

template<typename eigenArr , typename realParamT >

void mx::sigproc::frequencyGrid	(	eigenArr &	arr,
		realParamT	drT,
		eigenArr *	k_x,
		eigenArr *	k_y
	)

Create a 2-D frequency grid.

Definition at line 313 of file psdUtils.hpp.

References mx::sigproc::freq_sampling().

◆ frequencyGrid() [2/4]

template<typename eigenArr >

void mx::sigproc::frequencyGrid	(	eigenArr &	arr,
		typename eigenArr::Scalar	dt
	)

Create a frequency grid.

Definition at line 389 of file psdUtils.hpp.

References mx::sigproc::frequencyGrid().

◆ frequencyGrid() [3/4]

template<typename eigenArr >

void mx::sigproc::frequencyGrid	(	eigenArr &	arr,
		typename eigenArr::Scalar	dt,
		eigenArr &	k_x,
		eigenArr &	k_y
	)

Create a frequency grid.

Definition at line 396 of file psdUtils.hpp.

References mx::sigproc::frequencyGrid().

◆ frequencyGrid() [4/4]

template<typename realT , typename realParamT >

int mx::sigproc::frequencyGrid	(	std::vector< realT > &	vec,
		realParamT	dt,
		bool	fftOrder = `true`
	)

Create a 1-D frequency grid.

Template Parameters

realT	a real floating point type
realParamT	a real floating point type, convenience to avoid double-float confusion.

Test:: Verify creation of a 1D frequency grid [test doc]

Parameters

[out]	vec	vec the pre-allocated vector, on return contains the frequency grid
[in]	dt	dt the temporal sampling of the time series
[in]	fftOrder	fftOrder [optional] if true the frequency grid is in FFT order

Definition at line 256 of file psdUtils.hpp.

Referenced by mx::sigproc::frequencyGrid(), mx::sigproc::frequencyGrid(), SCENARIO(), and SCENARIO().

◆ kneePSD()

template<typename floatT >

int mx::sigproc::kneePSD	(	std::vector< floatT > &	psd,
		std::vector< floatT > &	f,
		floatT	beta,
		floatT	fn,
		floatT	alpha
	)

Generate a 1-D "knee" PSD.

Populates an Eigen array with

\[ P(f) = \frac{\beta}{ 1 + (f/f_n)^{\alpha}} \]

If you set \( T_0 \le 0 \) and \( t_0 = 0\) this reverts to a simple \( 1/f^\alpha \) law (i.e. it treats this as infinite outer scale and inner scale).

Template Parameters

floatT a floating point

Parameters

[out]	psd	the PSD vector, will be resized.
[in]	f	the frequency vector
[in]	beta	the scaling constant
[in]	fn	the knee frequency
[in]	alpha	the exponent, by convention \( alpha > 0 \).

Definition at line 699 of file psdUtils.hpp.

◆ normPSD() [1/2]

template<typename floatT , typename floatParamT >

floatT mx::sigproc::normPSD	(	Eigen::Array< floatT, Eigen::Dynamic, Eigen::Dynamic > &	psd,
		Eigen::Array< floatT, Eigen::Dynamic, Eigen::Dynamic > &	k,
		floatParamT	normT,
		floatT	kmin = `std::numeric_limits<floatT>::min()`,
		floatT	kmax = `std::numeric_limits<floatT>::max()`
	)

Normalize a 2-D PSD to have a given variance.

A frequency range can be specified for calculating the norm, otherwise the entire PSD is used. The entire PSD is normalized regardless.

Template Parameters

floatT	the floating point type of the PSD.
floatParamT	a floating point type, convenience to avoid double-float cofusion.

Test:: Verify scaling and normalization of augment1SidedPSD [test doc]

Parameters

	psd	[in.out] the PSD to normalize, will be altered.
[in]	k	the frequency grid for psd.
[in]	normT	the desired total variance (or integral) of the PSD.
[in]	kmin	[optiona] the minimum frequency of the range over which to normalize.
[in]	kmax	[optiona] the maximum frequency of the range over which to normalize.

Definition at line 497 of file psdUtils.hpp.

◆ normPSD() [2/2]

template<typename floatT , typename floatParamT >

int mx::sigproc::normPSD	(	std::vector< floatT > &	psd,
		std::vector< floatT > &	f,
		floatParamT	normT,
		floatT	fmin = `std::numeric_limits<floatT>::min()`,
		floatT	fmax = `std::numeric_limits<floatT>::max()`
	)

Normalize a 1-D PSD to have a given variance.

A frequency range can be specified to calculate the norm, otherwise f[0] to f[f.size()-1] is the range. The entire PSD is normalized regardless.

Template Parameters

floatT	the floating point type of the PSD.
floatParamT	a floating point type, convenience to avoid double-float cofusion.

Test:: Verify scaling and normalization of augment1SidedPSD [test doc]

Parameters

	psd	[in.out] the PSD to normalize, will be altered.
[in]	f	the frequency points for the PSD
[in]	normT	the desired total variance (or integral) of the PSD.
[in]	fmin	[optiona] the minimum frequency of the range over which to normalize.
[in]	fmax	[optiona] the maximum frequency of the range over which to normalize.

Definition at line 447 of file psdUtils.hpp.

Referenced by mx::AO::analysis::fourierTemporalPSD< _realT, aosysT >::intensityPSD(), SCENARIO(), and SCENARIO().

◆ oneoverf_norm()

template<typename realT >

realT mx::sigproc::oneoverf_norm	(	realT	fmin,
		realT	fmax,
		realT	alpha
	)

Calculate the normalization for a 1-D \( 1/|f|^\alpha \) PSD.

Parameters

[in]	fmin	is the minimum non-zero absolute value of frequency
[in]	fmax	is the maximum absolute value of frequencey
[in]	alpha	is the power-law exponent, by convention \( \alpha > 0 \).

Returns: the normalization for a 2-sided power law PSD.

Template Parameters

realT is the real floating point type used for calculations.

Definition at line 412 of file psdUtils.hpp.

Referenced by mx::sigproc::oneoverf_psd(), and mx::sigproc::vonKarmanPSD().

◆ oneoverf_psd()

template<typename eigenArrp , typename eigenArrf >

void mx::sigproc::oneoverf_psd	(	eigenArrp &	psd,
		eigenArrf &	freq,
		typename eigenArrp::Scalar	alpha,
		typename eigenArrp::Scalar	beta = `-1`
	)

Generates a \( 1/|f|^\alpha \) power spectrum.

Populates an Eigen array with

\[ P(|f| = 0) = 0 \]

\[ P(|f| > 0) = \frac{\beta}{|f|^{\alpha}} \]

Parameters

[out]	psd	is the array to populate
[in]	freq	is a frequency grid, must be the same logical size as psd
[in]	alpha	is the power law exponent, by convention \( alpha > 0 \).
[in]	beta	[optional is a normalization constant to multiply the raw spectrum by. If beta==-1 (default) then the PSD is normalized using oneoverf_norm.

Template Parameters

eigenArrp	is the Eigen-like array type of the psd
eigenArrf	is the Eigen-like array type of the frequency grid

Definition at line 580 of file psdUtils.hpp.

References mx::sigproc::oneoverf_norm().

◆ oneoverk_norm()

template<typename realT >

realT mx::sigproc::oneoverk_norm	(	realT	kmin,
		realT	kmax,
		realT	alpha
	)

Calculate the normalization for a 2-D \( 1/|k|^\alpha \) PSD.

Parameters

[in]	kmin	is the minimum non-zero absolute value of frequency
[in]	kmax	is the maximum absolute value of frequencey
[in]	alpha	is the power-law exponent, by convention \( \alpha > 0 \).

Returns: the normalization for a 2-D, 2-sided power law PSD.

Template Parameters

realT is the real floating point type used for calculations.

Definition at line 430 of file psdUtils.hpp.

◆ psdVar()

template<typename realT >

realT mx::sigproc::psdVar	(	const std::vector< realT > &	f,
		const std::vector< realT > &	PSD,
		realT	half = `0.5`
	)

Calculate the variance of a 1-D PSD.

By default uses trapezoid rule integration. This can be changed to mid-point integration.

If f.back() < 0, then a 2-sided PSD in FFT storage order is assumed. Otherwise, PSD is treated as 1-sided.

Returns: the variance of a PSD (the integral).

Template Parameters

realT the real floating point type

Test:: Scenario: calculating variance from a 1D PSD. [test doc]

Parameters

[in]	f	the frequency scale of the PSD.
[in]	PSD	the PSD to integrate.
[in]	half	[optional] controls if trapezoid (0.5) or mid-point (1.0) integration is used. Do not use other values.

Definition at line 136 of file psdUtils.hpp.

References mx::sigproc::psdVar1sided(), and mx::sigproc::psdVar2sided().

Referenced by mx::AO::analysis::fourierTemporalPSD< _realT, aosysT >::analyzePSDGrid(), mx::AO::analysis::fourierTemporalPSD< _realT, aosysT >::intensityPSD(), mx::sigproc::rebin1SidedPSD(), SCENARIO(), SCENARIO(), and mx::AO::analysis::speckleAmpPSD().

◆ psdVar1sided()

template<typename realT >

realT mx::sigproc::psdVar1sided	(	realT	df,
		const realT *	PSD,
		size_t	sz,
		realT	half = `0.5`
	)

Calculate the variance of a 1-D, 1-sided PSD.

By default uses trapezoid rule integration. This can be changed to mid-point integration.

Returns: the variance of a PSD (the integral).

Template Parameters

realT the real floating point type

Test:: Scenario: calculating variance from a 1D PSD [test doc]

Parameters

[in]	df	the frequency scale of the PSD
[in]	PSD	the PSD to integrate.
[in]	sz	the size of the PSD vector
[in]	half	[optional] controls if trapezoid (0.5) or mid-point (1.0) integration is used. Do not use other values.

Definition at line 66 of file psdUtils.hpp.

Referenced by mx::sigproc::averagePeriodogram< realT >::operator()(), and mx::sigproc::psdVar().

◆ psdVar2sided()

template<typename realT >

realT mx::sigproc::psdVar2sided	(	realT	df,
		const realT *	PSD,
		size_t	sz,
		realT	half = `0.5`
	)

Calculate the variance of a 1-D, 2-sided PSD.

By default uses trapezoid rule integration. This can be changed to mid-point integration.

Assumes the 2-sided PSD is in standard FFT storage order, and that sz is even.

Returns: the variance of a PSD (the integral).

Template Parameters

realT the real floating point type

Test:: Scenario: calculating variance from a 1D PSD. [test doc]

Parameters

[in]	df	the frequency scale of the PSD
[in]	PSD	the PSD to integrate.
[in]	sz	the size of the PSD vector
[in]	half	[optional] controls if trapezoid (0.5) or mid-point (1.0) integration is used. Do not use other values.

Definition at line 99 of file psdUtils.hpp.

Referenced by mx::sigproc::psdVar().

◆ psdVarDisabled()

template<typename eigenArrT >

eigenArrT::Scalar mx::sigproc::psdVarDisabled	(	eigenArrT &	freq,
		eigenArrT &	PSD,
		bool	trap = `true`
	)

Calculate the variance of a PSD.

By default uses trapezoid rule integration. This can be changed to mid-point integration.

This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

Returns: the variance of a PSD (the integral).

Template Parameters

realT the real floating point type

Parameters

[in]	freq	the frequency scale of the PSD
[in]	PSD	the PSD to integrate.
[in]	trap	[optional] controls if trapezoid (true) or mid-point (false) integration is used.

Definition at line 158 of file psdUtils.hpp.

◆ rebin1SidedPSD()

template<typename realT >

int mx::sigproc::rebin1SidedPSD	(	std::vector< realT > &	binFreq,
		std::vector< realT > &	binPSD,
		std::vector< realT > &	freq,
		std::vector< realT > &	PSD,
		realT	binSize,
		bool	binAtZero = `true`
	)

Rebin a PSD, including its frequency scale, to a larger frequency bin size (fewer bins)

The rebinning uses trapezoid integration within bins to ensure minimum signal loss.

Maintains DFT sampling. That is, if initial frequency grid is 0,0.1,0.2... and the binSize is 1.0, the new grid will be 0,1,2 (as opposed to 0.5, 1.5, 2.5).

This introduces a question of what to do with first half-bin, which includes 0. It can be integrated (binAtZero = true, the default). This may cause inaccurate behavior if the value of the PSD when f=0 is important (e.g. when analyzing correlated noise), so setting binAtZero=false causes the f=0 value to be copied (using the nearest neighbor if no f=0 point is in the input.

The last half bin is always integrated.

The output is variance normalized to match the input variance.

Template Parameters

realT the real floating point type

Parameters

[out]	binFreq	the binned frequency scale, resized.
[out]	binPSD	the binned PSD, resized.
[in]	freq	the frequency scale of the PSD to bin.
[in]	PSD	the PSD to bin.
[in]	binSize	in same units as freq
[in]	binAtZero	[optional] controls whether the zero point is binned or copied.

Definition at line 933 of file psdUtils.hpp.

References mx::sigproc::psdVar().

◆ vonKarmanPSD() [1/2]

template<typename eigenArrp , typename eigenArrf , typename alphaT , typename L0T , typename l0T , typename betaT >

void mx::sigproc::vonKarmanPSD	(	eigenArrp &	psd,
		eigenArrf &	freq,
		alphaT	alpha,
		L0T	L0 = `0`,
		l0T	l0 = `0`,
		betaT	beta = `-1`
	)

Generates a von Karman power spectrum.

Populates an Eigen array with

\[ P(k) = \frac{\beta}{ (k^2 + (1/L_0)^2)^{\alpha/2}} e^{ - k^2 l_0^2} \]

If you set \( L_0 \le 0 \) and \( l_0 = 0\) this reverts to a simple \( 1/f^\alpha \) law (i.e. it treats this as infinite outer scale and inner scale).

Parameters

[out]	psd	is the array to populate, allocated.
[in]	freq	is a frequency grid, must be the same logical size as psd
[in]	alpha	is the power law exponent, by convention \( alpha > 0 \).
[in]	L0	[optional] is the outer scale.
[in]	l0	[optional] is the inner scale.
[in]	beta	[optional] is a normalization constant to multiply the raw spectrum by. If beta==-1 (default) then the PSD is normalized using oneoverf_norm.

Template Parameters

eigenArrp	is the Eigen array type of the psd
eigenArrf	is the Eigen array type of the frequency grid

Definition at line 741 of file psdUtils.hpp.

References mx::sigproc::oneoverf_norm().

◆ vonKarmanPSD() [2/2]

template<typename floatT , typename floatfT , typename alphaT , typename T0T = double, typename t0T = double, typename betaT = double>

int mx::sigproc::vonKarmanPSD	(	std::vector< floatT > &	psd,
		std::vector< floatfT > &	f,
		alphaT	alpha,
		T0T	T0 = `0`,
		t0T	t0 = `0`,
		betaT	beta = `1`
	)

Generate a 1-D von Karman power spectrum.

Populates an Eigen array with

\[ P(f) = \frac{\beta}{ (f^2 + (1/T_0)^2)^{\alpha/2}} e^{ - f^2 t_0^2} \]

If you set \( T_0 \le 0 \) and \( t_0 = 0\) this reverts to a simple \( 1/f^\alpha \) law (i.e. it treats this as infinite outer scale and inner scale).

Template Parameters

floatT a floating point

Parameters

[out]	psd	the PSD vector, will be resized.
[in]	f	the frequency vector
[in]	alpha	the exponent, by convention \( alpha > 0 \).
[in]	T0	the outer scale, default is 0 (not used).
[in]	t0	the inner scale, default is 0 (not used).
[in]	beta	the scaling constant, default is 1

Definition at line 643 of file psdUtils.hpp.

Classes

Modules

Function Documentation

◆ augment1SidedPSD()

◆ augment1SidedPSDFreq()

◆ freq_sampling()

◆ frequency_grid1D()

◆ frequencyGrid() [1/4]

◆ frequencyGrid() [2/4]

◆ frequencyGrid() [3/4]

◆ frequencyGrid() [4/4]

◆ kneePSD()

◆ normPSD() [1/2]

◆ normPSD() [2/2]

◆ oneoverf_norm()

◆ oneoverf_psd()

◆ oneoverk_norm()

◆ psdVar()

◆ psdVar1sided()

◆ psdVar2sided()

◆ psdVarDisabled()

◆ rebin1SidedPSD()

◆ vonKarmanPSD() [1/2]

◆ vonKarmanPSD() [2/2]