Compute data spectra. More...
#include <Spectrum.h>
Public Member Functions | |
bool | AddData (const int aDataSize, const double *aData, const int aDataStart=0) |
Loads a data vector and updates the current PSD. | |
TGraph * | GetASD (const double aFrequencyMin, const double aFrequencyMax) |
Returns the current ASD as a TGraph. | |
TGraph * | GetASD (void) |
Returns the current ASD as a TGraph. | |
int | GetDataBufferLength (void) |
Returns the effective data buffer length [s] used to compute the PSD. | |
int | GetNSubSegmentsMax (const int aParity) |
Returns the maximum number of subsegments in the circular buffer. | |
double | GetPower (const double aFrequency) |
Returns the current power value at a given frequency. | |
TGraph * | GetPSD (const double aFrequencyMin, const double aFrequencyMax) |
Returns the current PSD as a TGraph between 2 frequencies. | |
TGraph * | GetPSD (void) |
Returns the current PSD as a TGraph. | |
double | GetSpectrumFrequency (const int aIndex) |
Returns the frequency of a given spectrum frequency index. | |
int | GetSpectrumNyquist (void) |
Returns the Nyquist frequency of the spectrum [Hz]. | |
double | GetSpectrumResolution (void) |
Returns the frequency resolution of the spectrum [Hz]. | |
int | GetSpectrumSize (void) |
Returns the number of points in the spectrum. | |
bool | GetStatus (void) |
Returns class status. | |
TGraph * | GetSubASD (const int aParity, const int aIndex) |
Returns the ASD for a given subsegment. | |
TGraph * | GetSubPSD (const int aParity, const int aIndex) |
Returns the PSD for a given subsegment. | |
bool | IsBufferEmpty (void) |
Checks the current status of the PSD buffer. | |
bool | LoadData (const int aDataSize, const double *aData, const int aDataStart=0) |
Loads a data vector and computes a new PSD. | |
void | Reset (void) |
Resets the current PSD and PSD circular buffer. | |
bool | WriteSubASD (const string aOutFileName) |
Writes all sub-segment ASDs in a ROOT file. | |
bool | WriteSubPSD (const string aOutFileName) |
Writes all sub-segment PSDs in a ROOT file. | |
Constructors and destructors | |
Spectrum (const int aSpectrumSize, const int aDataDuration, const int aDataSamplingFrequency, const int aVerbosity=0) | |
Constructor of the Spectrum class. | |
virtual | ~Spectrum (void) |
Destructor of the Spectrum class. |
Compute data spectra.
This class was designed to compute the one-sided power spectrum 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 cartoon below:
|--------------------------------| input data vector (time domain) |---|---|---|---|---|---|---|---| Neven 'even' subsegments |---|---|---|---|---|---|---| Nodd 'odd' subsegments
The input data vector is divided into Neven 'even' subsegments and Nodd 'odd' subsegments. Odd and even subsegments overlap by 50%. Neven and Nodd are determined by the size of the input vector and the requested frequency resolution for the PSD. The minimum requirement is to have at least Neven=1 and Nodd=0 to estimate the PSD (not the best estimate!). Each subsegment is hann-windowed and FFTed to get Nodd+Neven power spectra. The median PSDs are computed for odd and even subsegments separately and a median-to-average correction factor is applied. The 2 PSDs are finally averaged into one.
In fact, a dynamical PSD estimation is in use. This means that number of even/odd sub-segments is defined and fixed by the constructor of this class. Then, the PSD can be updated by providing an arbitrary number of subsegment data: see AddData().
The PSD is estimated with a given frequency resolution. A linear interpolation is used to compute the power for any given frequency.
Spectrum::Spectrum | ( | const int | aSpectrumSize, | |
const int | aDataDuration, | |||
const int | aDataSamplingFrequency, | |||
const int | aVerbosity = 0 | |||
) |
Constructor of the Spectrum class.
It defines the vector sizes and initiates the FFT plan (FFTW_MEASURE). The spectrum resolution is defined by the number of points used to compute the data spectrum (between 0 and the Nyquist frequency). A nominal duration over which the PSD is estimated must be given (see the method presented in the class description). From this a circular buffer of sub-segments is defined and will be used to apply the median-mean method.
The data sampling frequency and the spectrum size must be even numbers. The nominal data duration must be long enough to contain at least 1 subsegment (1 even and 0 odd), in other words, aDataDuration*aDataSamplingFrequency >= 2*aSpectrumSize.
The actual spectrum size is 'aSpectrumSize' + 1 to account for the DC frequency.
aSpectrumSize | number of points in the spectrum | |
aDataDuration | nominal data duration to estimate the PSD [s] | |
aDataSamplingFrequency | input data sampling frequency [Hz] | |
averbose | verbosity level |
virtual Spectrum::~Spectrum | ( | void | ) | [virtual] |
Destructor of the Spectrum class.
bool Spectrum::AddData | ( | const int | aDataSize, | |
const double * | aData, | |||
const int | aDataStart = 0 | |||
) |
Loads a data vector and updates the current PSD.
The data vector 'aData' (time domain) is used to update the current PSD. The method presented in this class description is applied.
The user must specify the size of the input vector to cover. It must be long enough to be able to compute at least 1 PSD (>=2*aSpectrumSize).
If the size of the input vector is smaller than the size specified in the constructor, the circular buffer is updated with this new data set. The oldest PSDs in the buffer are erased. If the size of the input vector is larger than the size specified in the constructor, only the end of the input vector is considered (down to the size specified in the constructor). This case is computationnally inefficient as sub-segmenmt PSDs are computed for nothing.
Optionally, the user can provide a starting index. The data vector will only be used starting at this index. The aDataSize value should then account for a shifted start (do not mess this up as it will not be checked!).
aDataSize | input vector size (number of samples) starting at aDataStart | |
aData | input data vector (time domain) | |
aDataStart | first data index to consider for the PSD estimation (=0 if from start) |
TGraph* Spectrum::GetASD | ( | const double | aFrequencyMin, | |
const double | aFrequencyMax | |||
) |
Returns the current ASD as a TGraph.
The returned TGraph must be deleted by the user. NULL is return if this function fails.
aFrequencyMin | minimum frequency value [Hz] | |
aFrequencyMax | minimum frequency value [Hz] |
TGraph* Spectrum::GetASD | ( | void | ) |
Returns the current ASD as a TGraph.
The returned TGraph must be deleted by the user. NULL is return if this function fails.
int Spectrum::GetDataBufferLength | ( | void | ) | [inline] |
Returns the effective data buffer length [s] used to compute the PSD.
int Spectrum::GetNSubSegmentsMax | ( | const int | aParity | ) | [inline] |
Returns the maximum number of subsegments in the circular buffer.
aParity | 0 for even subsegments, 1 for odd subsegments |
double Spectrum::GetPower | ( | const double | aFrequency | ) |
Returns the current power value at a given frequency.
The power value, is computed by linearly interpolating the PSD if the requested frequency is inside the PSD frequency range. If it is outside this range, the PSD is linearly extrapolated using the last 2 points in the range. If the extrapolation leads to a negative value, the power is forced to the last-known value (edge).
If the requested frequency is negative or zero, the power DC value is returned
-1.0 is returned if this function fails.
If this function is called many times in a row, it is more optimal if the requested frequencies are called in a sorted order (increasing or decreasing).
aFrequency | frequency value [Hz] |
TGraph* Spectrum::GetPSD | ( | const double | aFrequencyMin, | |
const double | aFrequencyMax | |||
) |
Returns the current PSD as a TGraph between 2 frequencies.
The returned TGraph must be deleted by the user. NULL is return if this function fails.
aFrequencyMin | minimum frequency value [Hz] | |
aFrequencyMax | minimum frequency value [Hz] |
TGraph* Spectrum::GetPSD | ( | void | ) |
Returns the current PSD as a TGraph.
The returned TGraph must be deleted by the user. NULL is return if this function fails.
double Spectrum::GetSpectrumFrequency | ( | const int | aIndex | ) | [inline] |
Returns the frequency of a given spectrum frequency index.
aIndex | frequency index (0=DC) |
int Spectrum::GetSpectrumNyquist | ( | void | ) | [inline] |
Returns the Nyquist frequency of the spectrum [Hz].
double Spectrum::GetSpectrumResolution | ( | void | ) | [inline] |
Returns the frequency resolution of the spectrum [Hz].
int Spectrum::GetSpectrumSize | ( | void | ) | [inline] |
Returns the number of points in the spectrum.
The actual spectrum size is the returned value -1 to remove the DC frequency.
bool Spectrum::GetStatus | ( | void | ) | [inline] |
Returns class status.
TGraph* Spectrum::GetSubASD | ( | const int | aParity, | |
const int | aIndex | |||
) |
Returns the ASD for a given subsegment.
The returned TGraph must be deleted by the user. NULL is return if this function fails.
aParity | 0 for even subsegment, 1 for odd subsegment | |
aIndex | subsegment index |
TGraph* Spectrum::GetSubPSD | ( | const int | aParity, | |
const int | aIndex | |||
) |
Returns the PSD for a given subsegment.
The returned TGraph must be deleted by the user. NULL is return if this function fails.
aParity | 0 for even subsegments, 1 for odd subsegments | |
aIndex | subsegment index |
bool Spectrum::IsBufferEmpty | ( | void | ) | [inline] |
Checks the current status of the PSD buffer.
Returns true if the PSD buffer is empty, i.e. not a single PSD is available.
bool Spectrum::LoadData | ( | const int | aDataSize, | |
const double * | aData, | |||
const 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.
aDataSize | input vector size (number of samples) starting at aDataStart | |
aData | input data vector (time domain) | |
aDataStart | first data index to consider for the PSD estimation (=0 if from start) |
void Spectrum::Reset | ( | void | ) |
Resets the current PSD and PSD circular buffer.
The PSD buffer is flushed. The current PSD is set to 0.
bool Spectrum::WriteSubASD | ( | const string | aOutFileName | ) |
Writes all sub-segment ASDs in a ROOT file.
aOutFileName | ROOT output file name |
bool Spectrum::WriteSubPSD | ( | const string | aOutFileName | ) |
Writes all sub-segment PSDs in a ROOT file.
aOutFileName | ROOT output file name |