DFT

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

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.

Functions
DAS_API const double *	Dft_getReal (Das2Dft *pThis, size_t *pLen)
	Return the real component after a calculation. More...
DAS_API Das2Dft *	new_Dft (DftPlan *pPlan, const char *sWindow)
	Create a new DFT calculator. More...
DAS_API void	del_Dft (Das2Dft *pThis)
	Free a DFT (Discrete Fourier Transform) calculator. More...
DAS_API DasErrCode	Dft_calculate (Das2Dft *pThis, const double *pReal, const double *pImg)
	Calculate a discrete Fourier transform. More...
DAS_API const double *	Dft_getImg (Das2Dft *pThis, size_t *pLen)
	Return the imaginary component after a calculation. More...
DAS_API const double *	Dft_getMagnitude (Das2Dft *pThis, size_t *pLen)
	Get the amplitude magnitude vector from a calculation. More...
DAS_API Das2Psd *	new_Psd (DftPlan *pPlan, bool bCenter, const char *sWindow)
	Create a new Power Spectral Density Calculator. More...
DAS_API void	del_Das2Psd (Das2Psd *pThis)
	Free a Power Spectral Density calculator. More...
DAS_API DasErrCode	Psd_calculate (Das2Psd *pThis, const double *pReal, const double *pImg)
	Calculate a Power Spectral Density (periodogram) More...
DAS_API DasErrCode	Psd_calculate_f (Das2Psd *pThis, const float *pReal, const float *pImg)
	The floating point array input analog of Psd_calaculate() More...
DAS_API double	Psd_powerRatio (const Das2Psd *pThis, double *pInput, double *pOutput)
	Provide a comparison of the input power and the output power. More...
DAS_API const double *	Psd_get (const Das2Psd *pThis, size_t *pLen)
	Get the amplitude magnitude vector from a calculation. 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 ...
  
}

Function Documentation

Dft_getReal()

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

Return the real component after a calculation.

Parameters

pThis
pLen

Returns

new_Dft()

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

del_Dft()

DAS_API 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

Dft_calculate()

DAS_API 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

Dft_getImg()

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

Return the imaginary component after a calculation.

Parameters

pThis
pLen

Returns

Dft_getMagnitude()

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

new_Psd()

DAS_API 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

del_Das2Psd()

DAS_API 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

Psd_calculate()

DAS_API 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

Psd_calculate_f()

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

Psd_powerRatio()

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

Psd_get()

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