Compute the noise power spectral density. More...

#include <Spectrum.h>

Collaboration diagram for Spectrum:

Public Member Functions
bool	AddData (const unsigned int aDataSize, const double *aData, const unsigned int aDataStart=0)
	Loads a data vector and updates the current PSD. More...

double	GetAmplitude (const double aFrequency)
	Returns the current ASD value at a given frequency. More...

TGraph *	GetAmplitudePeriodogram (const unsigned int aParity, const unsigned int aIndex)
	Returns a copy of a given (amplitude) periodogram. More...

TGraph *	GetASD (const double aFrequencyMin, const double aFrequencyMax)
	Returns a copy of the current ASD as a TGraph between 2 frequencies. More...

TGraph *	GetASD (void)
	Returns a copy of the current ASD as a TGraph. More...

unsigned int	GetDataBufferLength (void)
	Returns the effective circular buffer length [s] used to compute the PSD. More...

unsigned int	GetNSubSegmentsMax (const unsigned int aParity)
	Returns the maximum number of sub-segments in the buffer. More...

TGraph *	GetPeriodogram (const unsigned int aParity, const unsigned int aIndex)
	Returns a copy of a given periodogram. More...

double	GetPower (const double aFrequency)
	Returns the current PSD value at a given frequency. More...

double	GetPower (const unsigned int aFrequencyIndex)
	Returns the current PSD value at a given frequency index. More...

TGraph *	GetPSD (const double aFrequencyMin, const double aFrequencyMax)
	Returns a copy of the current PSD as a TGraph between 2 frequencies. More...

TGraph *	GetPSD (void)
	Returns a copy of the current PSD as a TGraph. More...

double	GetSpectrumFrequency (const unsigned int aIndex)
	Returns the frequency of a given spectrum frequency index [Hz]. More...

unsigned int	GetSpectrumNyquist (void)
	Returns the Nyquist frequency of the spectrum [Hz]. More...

double	GetSpectrumResolution (void)
	Returns the frequency resolution of the spectrum [Hz]: \(\delta f\). More...

unsigned int	GetSpectrumSize (void)
	Returns the number of points in the spectrum: \(n\). More...

bool	LoadData (const unsigned int aDataSize, const double *aData, const unsigned int aDataStart=0)
	Loads a data vector and computes a new PSD. More...

void	Reset (void)
	Resets the current PSD and the periodogram circular buffer. More...

bool	WritePeriodogram (const string aOutFileName)
	Writes periodograms and the resulting PSD in a ROOT file. More...

Constructors and destructors
	Spectrum (const unsigned int aSpectrumSize, const unsigned int aDataDuration, const unsigned int aDataSamplingFrequency, const unsigned int aVerbosity=0)
	Constructor of the Spectrum class. More...

virtual	~Spectrum (void)
	Destructor of the Spectrum class. More...

Private Member Functions
void	ComputeMedianMean (void)
	Computes median-mean PSD. More...

Private Attributes
fft *	fft_psd
	FFT plan to compute the PSD. More...

unsigned int	fVerbosity
	Verbosity level. More...

double *	HannWindow
	Hann window vector. More...

unsigned int	nSubSegments [2]
	Number of even/odd sub-subsegments. More...

double **	periodogram [2]
	Buffer for periodograms. More...

unsigned int	sampling
	Data sampling frequency [Hz]. More...

unsigned int	write_index [2]
	Write index in the buffer. More...

Detailed Description

Compute the noise power spectral density.

This class is designed to estimate the one-sided noise power spectral density of a given stretch of data (PSD). It uses the median-mean method described in gr-qc/0509116. This method can be represented by the diagram below:

Data segments to estimate the noise PSD.

The input data is divided into \(N_e\) "even" subsegments and \(N_o\) "odd" subsegments (Spectrum::Spectrum()) of duration \(T_s\). Odd and even subsegments overlap by 50%. The number of subsegments, \(N_e\) and \(N_o\), is given by the duration of the input data \(T\) and the requested frequency resolution for the PSD. The input data is sampled at a frequency \(f_s\); therefore the total number of data points is \(N = T\times f_s\). The PSD frequency resolution \(\delta f\) is defined by the number of points, \(n = f_s/(2\delta f)\), to cover frequencies from 0 Hz to the Nyquist frequency, \(f_s/2\). We have:

\[ \begin{cases} N_e = \lfloor{ \frac{T}{T_s} }\rfloor = \lfloor{ \frac{N}{2n} }\rfloor = \lfloor{ \frac{N\delta f}{f_s} } \rfloor \\ N_o = \lfloor{ \frac{T-T_s/2}{T_s} }\rfloor = \lfloor{ \frac{N-n}{2n} }\rfloor \end{cases} \]

Note: Some input data may be left unused at the end of the data segment (in red on the diagram above).

Warning: The spectrum frequency resolution \(n\) and the input data length \(T\) must be chosen to have at least one subsegment.

Each subsegment \(d_i(t)\) is hann-windowed with \(w(t)\) and Fourier transformed: \(\widetilde{(wd_i)}(f)\). Then, periodograms \(P_i(f)\) are evaluated for each subsegment:

\[ P_i(f) = \frac{2}{T_s}\left| \widetilde{(wd_i)}(f)\right|^2 \]

A median PSD is computed using periodograms for odd and even subsegments separately and a median-to-mean correction factor is applied (medianbiasfactor()). The 2 PSDs are finally averaged into one.

In fact, a dynamical PSD estimation is implemented here. The number of even and odd sub-segments is defined and fixed in the constructor of this class (Spectrum::Spectrum()). This defines a circular buffer of periodograms. Then, segments of data are provided with the AddData() method. This segment of data is used to compute periodograms, fill the circular buffer, and update the PSD estimate.

Author: Florent Robinet

Constructor & Destructor Documentation

◆ Spectrum()

Spectrum::Spectrum	(	const unsigned int	aSpectrumSize,
		const unsigned int	aDataDuration,
		const unsigned int	aDataSamplingFrequency,
		const unsigned int	aVerbosity = `0`
	)

Constructor of the Spectrum class.

It defines the data sizes and initiates the FFT plan (FFTW_MEASURE). The spectrum resolution is defined by the number of points used to compute the data spectrum: \(n\). A nominal duration over which the PSD is estimated must be given: \(T\). A circular buffer of sub-segments is defined and is used to apply the median-mean method.

Parameters

[in]	aSpectrumSize	Number of points in the spectrum (one sided): \(n\).
[in]	aDataDuration	Nominal data duration to estimate the PSD [s]: \(T\).
[in]	aDataSamplingFrequency	Input data sampling frequency [Hz]: \(f_s\).
[in]	aVerbosity	Verbosity level.

◆ ~Spectrum()

Spectrum::~Spectrum ( void )

virtual

Destructor of the Spectrum class.

Member Function Documentation

◆ AddData()

bool Spectrum::AddData	(	const unsigned int	aDataSize,
		const double *	aData,
		const unsigned int	aDataStart = `0`
	)

Loads a data vector and updates the current PSD.

The data 'aData' is used to update the current PSD. The user must specify the size of the input data. If the size of the input data is smaller than the size specified in the constructor, the circular buffer is updated with this new data set. The oldest periodograms in the buffer are removed. If the size of the input data is larger than the size specified in the constructor, only the end of the input data is considered.

Note: Optionally, the user can provide a start index. The data vector is scanned from this index.

The input data is broken into odd and even subsegments as described in this class introduction. The periodograms are computed and saved in the circular buffer. Finally, a new PSD is estimated.

Returns: true if the input data is correctly used for the PSD estimate. false if the input data is of too short duration to compute the PSD.

Parameters

[in]	aDataSize	Input vector size (number of points) starting at aDataStart.
[in]	aData	Input time series.
[in]	aDataStart	First data index to consider for the PSD estimation.

◆ ComputeMedianMean()

void Spectrum::ComputeMedianMean ( void )

private

Computes median-mean PSD.

◆ GetAmplitude()

double Spectrum::GetAmplitude ( const double aFrequency )

inline

Returns the current ASD value at a given frequency.

The ASD value is the square root of the PSD returned with GetPower().

Parameters

[in] aFrequency Frequency value [Hz].

◆ GetAmplitudePeriodogram()

TGraph * Spectrum::GetAmplitudePeriodogram	(	const unsigned int	aParity,
		const unsigned int	aIndex
	)

Returns a copy of a given (amplitude) periodogram.

Warning: The returned TGraph must be deleted by the user.

Parameters

[in]	aParity	0 for even subsegments, 1 for odd subsegments.
[in]	aIndex	Subsegment index.

◆ GetASD() [1/2]

TGraph * Spectrum::GetASD	(	const double	aFrequencyMin,
		const double	aFrequencyMax
	)

Returns a copy of the current ASD as a TGraph between 2 frequencies.

Warning: The returned TGraph must be deleted by the user.

Parameters

[in]	aFrequencyMin	Minimum frequency value [Hz].
[in]	aFrequencyMax	Maximum frequency value [Hz].

◆ GetASD() [2/2]

TGraph * Spectrum::GetASD ( void )

Returns a copy of the current ASD as a TGraph.

Warning: The returned TGraph must be deleted by the user.

◆ GetDataBufferLength()

unsigned int Spectrum::GetDataBufferLength ( void )

inline

Returns the effective circular buffer length [s] used to compute the PSD.

◆ GetNSubSegmentsMax()

unsigned int Spectrum::GetNSubSegmentsMax ( const unsigned int aParity )

inline

Returns the maximum number of sub-segments in the buffer.

Parameters

[in] aParity 0 for even subsegments, 1 for odd subsegments

◆ GetPeriodogram()

TGraph * Spectrum::GetPeriodogram	(	const unsigned int	aParity,
		const unsigned int	aIndex
	)

Returns a copy of a given periodogram.

Warning: The returned TGraph must be deleted by the user.

Parameters

[in]	aParity	0 for even subsegments, 1 for odd subsegments.
[in]	aIndex	Subsegment index.

◆ GetPower() [1/2]

double Spectrum::GetPower ( const double aFrequency )

Returns the current PSD value at a given frequency.

The power value is computed by interpolating the PSD if the requested frequency is inside the PSD frequency range (0 to Nyquist). If it is outside this range, the PSD is extrapolated.

Note: If the requested frequency is negative or zero, the DC power value is returned.

Parameters

[in] aFrequency Frequency value [Hz].

◆ GetPower() [2/2]

double Spectrum::GetPower ( const unsigned int aFrequencyIndex )

inline

Returns the current PSD value at a given frequency index.

Parameters

aFrequencyIndex Frequency index.

Precondition: The frequency index must be a valid index: from 0 to the value returned by GetSpectrumSize().

◆ GetPSD() [1/2]

TGraph * Spectrum::GetPSD	(	const double	aFrequencyMin,
		const double	aFrequencyMax
	)

Returns a copy of the current PSD as a TGraph between 2 frequencies.

Warning: The returned TGraph must be deleted by the user.

Parameters

[in]	aFrequencyMin	Minimum frequency value [Hz].
[in]	aFrequencyMax	Maximum frequency value [Hz].

◆ GetPSD() [2/2]

TGraph * Spectrum::GetPSD ( void )

Returns a copy of the current PSD as a TGraph.

Warning: The returned TGraph must be deleted by the user.

◆ GetSpectrumFrequency()

double Spectrum::GetSpectrumFrequency ( const unsigned int aIndex )

inline

Returns the frequency of a given spectrum frequency index [Hz].

Parameters

[in] aIndex Frequency index (0=DC).

◆ GetSpectrumNyquist()

unsigned int Spectrum::GetSpectrumNyquist ( void )

inline

Returns the Nyquist frequency of the spectrum [Hz].

◆ GetSpectrumResolution()

double Spectrum::GetSpectrumResolution ( void )

inline

Returns the frequency resolution of the spectrum [Hz]: \(\delta f\).

◆ GetSpectrumSize()

unsigned int Spectrum::GetSpectrumSize ( void )

inline

Returns the number of points in the spectrum: \(n\).

◆ LoadData()

bool Spectrum::LoadData	(	const unsigned int	aDataSize,
		const double *	aData,
		const unsigned int	aDataStart = `0`
	)

Loads a data vector and computes a new PSD.

This function is the same as AddData() but the PSD circular buffer is flushed before loading the new data vector.

Returns: true if the input data is correctly used for the PSD estimate. false if the input data is of too short duration to compute the PSD.

Parameters

[in]	aDataSize	Input vector size (number of samples) starting at aDataStart
[in]	aData	Input time series.
[in]	aDataStart	First data index to consider for the PSD estimation.

◆ Reset()

void Spectrum::Reset ( void )

Resets the current PSD and the periodogram circular buffer.

The periodogram buffer is flushed. The current PSD is set to 0.

◆ WritePeriodogram()

bool Spectrum::WritePeriodogram ( const string aOutFileName )

Writes periodograms and the resulting PSD in a ROOT file.

Parameters

[in] aOutFileName ROOT output file name (use a ".root" extension).

Returns: false if the file can not be opened, and true otherwise.

Member Data Documentation

◆ fft_psd

fft* Spectrum::fft_psd

private

FFT plan to compute the PSD.

◆ fVerbosity

unsigned int Spectrum::fVerbosity

private

Verbosity level.

◆ HannWindow

double* Spectrum::HannWindow

private

Hann window vector.

◆ nSubSegments

unsigned int Spectrum::nSubSegments[2]

private

Number of even/odd sub-subsegments.

◆ periodogram

double** Spectrum::periodogram[2]

private

Buffer for periodograms.

◆ sampling

unsigned int Spectrum::sampling

private

Data sampling frequency [Hz].

◆ write_index

unsigned int Spectrum::write_index[2]

private

Write index in the buffer.

The documentation for this class was generated from the following files:

src/Streams/Spectrum.h
src/Streams/Spectrum.cc

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Spectrum()

◆ ~Spectrum()

Member Function Documentation

◆ AddData()

◆ ComputeMedianMean()

◆ GetAmplitude()

◆ GetAmplitudePeriodogram()

◆ GetASD() [1/2]

◆ GetASD() [2/2]

◆ GetDataBufferLength()

◆ GetNSubSegmentsMax()

◆ GetPeriodogram()

◆ GetPower() [1/2]

◆ GetPower() [2/2]

◆ GetPSD() [1/2]

◆ GetPSD() [2/2]

◆ GetSpectrumFrequency()

◆ GetSpectrumNyquist()

◆ GetSpectrumResolution()

◆ GetSpectrumSize()

◆ LoadData()

◆ Reset()

◆ WritePeriodogram()

Member Data Documentation

◆ fft_psd

◆ fVerbosity

◆ HannWindow

◆ nSubSegments

◆ periodogram

◆ sampling

◆ write_index