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...

Topics
	Autoregressive Power Spectra
	PSD Filter
	Filtering with a PSD to generate correlated noise.
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.
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.
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.
template<typename eigenArrT>
eigenArrT::Scalar	mx::sigproc::psdVarDisabled (eigenArrT &freq, eigenArrT &PSD, bool trap=true)
	Calculate the variance of a PSD.
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.
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<typename eigenArr, typename realParamT>
void	mx::sigproc::frequencyGrid (eigenArr &arr, realParamT drT, eigenArr k_x, eigenArr k_y)
	Create a 2-D frequency grid.
template<typename eigenArr>
void	mx::sigproc::frequencyGrid (eigenArr &arr, typename eigenArr::Scalar dt)
	Create a frequency grid.
template<typename eigenArr>
void	mx::sigproc::frequencyGrid (eigenArr &arr, typename eigenArr::Scalar dt, eigenArr &k_x, eigenArr &k_y)
	Create a frequency grid.
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.
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.
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.
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.
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.
template<typename floatT, typename floatfT, typename alphaT, typename T0T = double, typename t0T = double, typename betaT = double>
mx::error_t	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.
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.
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.
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.
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.
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).

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.

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 827 of file psdUtils.hpp.

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

◆ 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}

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 885 of file psdUtils.hpp.

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

◆ 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 199 of file psdUtils.hpp.

Referenced by frequencyGrid().

◆ 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 315 of file psdUtils.hpp.

References 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 391 of file psdUtils.hpp.

References 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 398 of file psdUtils.hpp.

References 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.

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 258 of file psdUtils.hpp.

References mx::invalidarg, and mx::internal::mxlib_error_report().

Referenced by frequencyGrid(), frequencyGrid(), TEST_CASE(), and TEST_CASE().

◆ 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 710 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.

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.

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 448 of file psdUtils.hpp.

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

◆ 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 414 of file psdUtils.hpp.

Referenced by oneoverf_psd(), and 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 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 432 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

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 135 of file psdUtils.hpp.

References psdVar1sided(), and psdVar2sided().

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

◆ 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

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 65 of file psdUtils.hpp.

Referenced by mx::sigproc::averagePeriodogram< realT >::operator()(), and 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

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 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 161 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 942 of file psdUtils.hpp.

References 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 752 of file psdUtils.hpp.

References oneoverf_norm().

◆ vonKarmanPSD() [2/2]

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

mx::error_t 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).

Returns: error_t::noerror on success

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 644 of file psdUtils.hpp.

References mx::noerror.

Classes

Topics

Function Documentation

◆ augment1SidedPSD()

◆ augment1SidedPSDFreq()

◆ freq_sampling()

◆ 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]