Discrete Fourier transforms and power spectral density estimation. More...

Detailed Description

Discrete Fourier transforms and power spectral density estimation.

This module provides an amplitude preserving 1-D Fourier transform and power preserving Power Spectral Density estimator. One key concept for this module is FFT plans can not be created and destroyed while transforms are running. On Linux the FFTW3 library is used to provide fast transform capability, the library to use for Windows is yet to be determined.

On linux if the file /etc/fftw/wisdom exists it will be loaded during the call to das_init(). Pre-planning FFT operations can significantly increase the speed of new_DftPlan() calls

The following example uses the pthread library to demonstrate running four simultaneous transforms.

* 
* #define SEC_IN_DAY 86400
* #define FREQUENCIES (SEC_IN_DAY/2 + 1)
* 
* // I/O structure for worker threads
* struct thread_data {
*   DftPlan* pPlan;
*   double* pInput;
*   double* pOutput;
* };
* 
* void* doTransform(void* vpThreadData){
*   struct thread_data* pTd = (struct thread_data*)vpThreadData;
* 
*   // DFT objects should be created on a per-thread basis
*   Das2Dft* pDft = new_Dft(pTd->pPlan, "HANN");
*   
*   // Run the transform saving the complex output to DFT obj memory
*   Dft_calculate(pDft, pTd->pInput, NULL);
*
*   // Copy magnitude out to the location designated be the main thread
*   size_t* uLen = 0;
*   const double* pMag = Dft_getMagnitude(pDft, &uLen);
*   assert(uLen == FREQUENCIES);
*   memcpy(pTd->pOutput, pMag, uLen*sizeof(double));
* 
*   del_Dft(pDft);
*   return NULL;
* }
* 
* main(){
*   
*   // Initialize libdas2
*    das_init(0, 0, DAS_LL_NOTICE, NULL);
* 
*   double* timeseries = (double*) calloc(SEC_IN_DAY*4, sizeof(double));
*   
*   // Do something to fill time array ...
*  
*   double* spectra = (double*) calloc(FREQUENCIES*4, sizeof(double));
* 
*   // Make a plan for calculating the DFTs
*   DftPlan* pPlan = new_DftPlan(SEC_IN_DAY, DAS2DFT_FORWARD, "my_program");
* 
*   // Now calculate 4 spectra at the same time
*   struct thread_data[4] = {NULL};
*   pthread_t threads[4];
*   for(int i = 0; i < 4; ++i){
*      thread_data[i].pPlan = pPlan;
*      thread_data[i].pInput = timeseries + SEC_IN_DAY*i;
*      thread_data[i].pOutput = spectra + FREQUENCIES*i;
*      pthread_create(threads+i, NULL, doTransform, thread_data+i);
*   }
*   
*   for(int i = 0; i < 4; ++i) pthread_join(threads + i);
*   
*   del_DftPlan(pPlan);
*   
*   // Do something with all 4 spectra ...
*   
* }
* 
* 

Data Structures
struct	Das2Dft
	An amplitude preserving Discrete Fourier Transform converter. More...

struct	Das2Psd
	A power spectral density estimator (periodogram) More...

Typedefs
typedef struct dft_plan	DftPlan
	An structure containing a set of global planning data for DFTs to be preformed. More...

Functions
const double *	Dft_getReal (Das2Dft pThis, size_t pLen)
	Return the real component after a calculation. More...

Das2Dft *	new_Dft (DftPlan pPlan, const char sWindow)
	Create a new DFT calculator. More...

void	del_Dft (Das2Dft *pThis)
	Free a DFT (Discrete Fourier Transform) calculator. More...

DasErrCode	Dft_calculate (Das2Dft pThis, const double pReal, const double *pImg)
	Calculate a discrete Fourier transform. More...

const double *	Dft_getImg (Das2Dft pThis, size_t pLen)
	Return the imaginary component after a calculation. More...

const double *	Dft_getMagnitude (Das2Dft pThis, size_t pLen)
	Get the amplitude magnitude vector from a calculation. More...

Das2Psd *	new_Psd (DftPlan pPlan, bool bCenter, const char sWindow)
	Create a new Power Spectral Density Calculator. More...

void	del_Das2Psd (Das2Psd *pThis)
	Free a Power Spectral Density calculator. More...

DasErrCode	Psd_calculate (Das2Psd pThis, const double pReal, const double *pImg)
	Calculate a Power Spectral Density (periodogram) More...

DasErrCode	Psd_calculate_f (Das2Psd pThis, const float pReal, const float *pImg)
	The floating point array input analog of Psd_calaculate() More...

double	Psd_powerRatio (const Das2Psd pThis, double pInput, double *pOutput)
	Provide a comparison of the input power and the output power. More...

const double *	Psd_get (const Das2Psd pThis, size_t pLen)
	Get the amplitude magnitude vector from a calculation. More...

Typedef Documentation

typedef struct dft_plan DftPlan

An structure containing a set of global planning data for DFTs to be preformed.

Function Documentation

const double* Dft_getReal	(	Das2Dft *	pThis,
		size_t *	pLen
	)

Return the real component after a calculation.

Parameters

pThis
pLen

Returns

Das2Dft * new_Dft	(	DftPlan *	pPlan,
		const char *	sWindow
	)

Create a new DFT calculator.

This function allocates re-usable storage for Fourier transform output, but in order to preform calculations a re-usable plan object must be provided

Parameters

pPlan	- A Transform plan object. The reference count of DFTs using this plan will be incremented.
sWindow	- A named window to apply to the data. If NULL then no window will be used.

Returns: A new Das2Dft object allocated on the heap.

void del_Dft ( Das2Dft * pThis )

Free a DFT (Discrete Fourier Transform) calculator.

Parameters

pThis the DFT calculator to free, the caller should set the object pointer to NULL after this call. Calling this also deletes the reference count for the associated DftPlan object

DasErrCode Dft_calculate	(	Das2Dft *	pThis,
		const double *	pReal,
		const double *	pImg
	)

Calculate a discrete Fourier transform.

Using the calculation plan setup in the constructor, calculate a discrete Fourier transform. When this function is called internal storage of any previous DFT calculations (if any) are over written.

Parameters

pThis	The DFT object which will hold the result memory
pReal	An input vector of with the same length as the plan object provided to the constructor
pImg	The imaginary (or quadrature phase) input vector with the same length as pRead. For a purely real signal this vector is NULL.

Returns: 0 (DAS_OKAY) if the calculation was successful, a non-zero error code otherwise

const double * Dft_getImg	(	Das2Dft *	pThis,
		size_t *	pLen
	)

Return the imaginary component after a calculation.

Parameters

pThis
pLen

Returns

const double * Dft_getMagnitude	(	Das2Dft *	pThis,
		size_t *	pLen
	)

Get the amplitude magnitude vector from a calculation.

Scale the stored DFT so that it preserves amplitude, and get the magnitude. For real-valued inputs (complex pointer = 0) the 'positive' and 'negative' frequencies are combined. For complex input vectors this is not the case since all DFT output amplitudes are unique. Stated another way, for complex input signals components above the Nyquist frequency have meaningful information.

Parameters

pThis	The DFT calculator object which has previously been called to calculate a result.
pLen	The vector length. In general this is NOT the same as the input time series length. For real-value input signals (complex input is NULL, this is N/2 + 1. For complex input signals this is N.

Returns: A pointer to an internal holding bin for the real signal magnitude values.

Warning: If Dft_calculate() is called again, the return pointer can be invalidated. If a permanent result is needed after subsequent Dft_calculate() calls, copy these data to another buffer.

Das2Psd * new_Psd	(	DftPlan *	pPlan,
		bool	bCenter,
		const char *	sWindow
	)

Create a new Power Spectral Density Calculator.

This estimator uses the equations given in Numerical Recipes in C, section 13.4, but not any of the actual Numerical Recipes source code.

Parameters

pPlan	- A Transform plan object. The reference count of DFTs using this plan will be incremented.
bCenter	If true, input values will be centered on the Mean value. This shifts-out the DC component from the input.
sWindow	A named window to use for the data. Possible values are: "hann" - Use a hann window as defined at http://en.wikipedia.org/wiki/Hann_function NULL - Use a square window. (i.e. 'multiply' all data by 1.0)

Returns: A new Power Spectral Density estimator allocated on the heap

void del_Das2Psd ( Das2Psd * pThis )

Free a Power Spectral Density calculator.

Parameters

pThis the PSD calculator to free, the caller should set the object pointer to NULL after this call. Calling this also deletes the reference count for the associated DftPlan object

DasErrCode Psd_calculate	(	Das2Psd *	pThis,
		const double *	pReal,
		const double *	pImg
	)

Calculate a Power Spectral Density (periodogram)

Using the calculation plan setup in the constructor, calculate a discrete Fourier transform. When this function is called internal storage of any previous DFT calculations (if any) are over written.

Parameters

pThis	The PSD calculator object
pReal	An input vector of with the same length as the plan object provided to the constructor
pImg	The imaginary (or quadrature phase) input vector with the same length as pRead. For a purely real signal this vector is NULL.

Returns: 0 (DAS_OKAY) if the calculation was successful, a non-zero error code otherwise

DasErrCode Psd_calculate_f	(	Das2Psd *	pThis,
		const float *	pReal,
		const float *	pImg
	)

The floating point array input analog of Psd_calaculate()

Internal calculations are still handled in double precision.

double Psd_powerRatio	(	const Das2Psd *	pThis,
		double *	pInput,
		double *	pOutput
	)

Provide a comparison of the input power and the output power.

During the Psd_calculate() call the average magnitude of the input vector is saved along with the average magnitude of the output vector (divided by the Window summed and squared). These two measures of power should always be close to each other when using a hann window. When using a NULL window they should be almost identical, to within rounding error. The two measures are:

             N-1
         1  ----   2      2
 Pin =  --- \    r    +  i
         N  /     n       n
            ----
             n=0

               N-1
          1   ----   2      2
 Pout =  ---  \    R    +  I
         Wss  /     k       k
              ----
               k=0

where Wss collapses to N**2 when a NULL (square) window is used. The reason that the Pout has an extra factor of N in the denominator is due to the following identity for the discrete Fourier transform (Parseval's theorem):

    N-1                   N-1
   ----   2    2      1  ----  2    2
   \     r  + i   =  --- \    R  + I
   /      n    n      N  /     n    n
   ----                  ----
    n=0                   k=0

Where r and i are the real and imaginary input amplitudes, and R and I are the DFT real and imaginary output values.

Parameters

pThis	A PSD calculator for which Psd_calculate has been called
pInput	A pointer to store the input power. If NULL, the input power will no be saved separately.
pOutput	A pointer to store the output power. If NULL, the output power will no be saved separately.

Returns: The ratio of Power Out divided by Power In. (Gain).

const double * Psd_get	(	const Das2Psd *	pThis,
		size_t *	pLen
	)

Get the amplitude magnitude vector from a calculation.

Scale the stored DFT so that it preserves amplitude, and get the magnitude. For real-valued inputs (complex pointer = 0) the 'positive' and 'negative' frequencies are combined. For complex input vectors this is not the case since all DFT output amplitudes are unique. Stated another way, for complex input signals components above the Nyquist frequency have meaningful information.

Parameters

pThis	The DFT calculator object which has previously been called to calculate a result.
pLen	The vector length. In general this is NOT the same as the input time series length. For real-value input signals (complex input is NULL, this is N/2 + 1. For complex input signals this is N.

Returns: A pointer to an internal holding bin for the real signal magnitude values.

Warning: If Psd_calculate() is called again, the return pointer can be invalidated. If a permanent result is needed after subsequent Psd_calculate() calls, copy these data to another buffer.

Detailed Description

Data Structures

Typedefs

Functions

Typedef Documentation

Function Documentation