Function: avg_spec()

0 void avg_spec(float *data,float *average,int npoint,int *reset,float srate,float decaytime,int windowtype,int overlap)

This routine calculates the power spectrum of the (time-domain) input stream data[ ], averaged over time with a user-set exponential decay, several possible choices of windowing and the possibility to overlap data in subsequent calls.

The arguments are:

data: Input. The time domain input samples are contained in data[0..N-1], with the data sample at time $t= n \Delta t$ contained in data[n].

average: Output. The one sided power spectrum is returned in average[0,..N/2-1]. The value of average[m] is the average power spectrum at frequency

$$f={{\tt m} \times {\tt srate} \over {\tt N}}.$$ (16.3.320)

We do not output the value of the average at the Nyquist frequency, which would be the (non-existent) array element average[N]. The units of average[ ] are ${\tt data[ ]}^2/{\rm Hz}$. Note: the elements of average[ ] must not be changed in between successive calls to avg_spec().

npoint: Input. The number of points ${\tt npoint}=N$ input. This must be an integer power of two.

reset: Input. If set to zero, then any past contribution to the average power spectrum is initialized to zero, and a new average is begun with the current input data.

srate: Input. The sample rate $1/\Delta t$ of the input data, in Hz.

decaytime: Input. The characteristic (positive) decay time $\tau$ in seconds, to use for the moving (exponentially-decaying) average described below. If no averaging over time is wanted, simply set decaytime to be small compared to $N \Delta t$.

windowtype: Input. Sets the type of window used in power spectrum estimation. Rectangular windowing (i.e., no windowing) is windowtype=0, Hann windowing is windowtype=1, Welch windowing is windowtype=2 and Bartlett windowing is windowtype=3. See [1] for a discussion of windowing and the definitions of these window types.

overlap: Input. Must be either zero or unity. If set to unity, then the data is overlapped by $N/2$ points with previous data (see below for a description of the overlapping procedure). When set to zero no overlapping is performed.

The methods used in this routine are quite similar to those used in the Numerical Recipes [1] routine spctrm(), and the reader interested in the details of this routine should first read the corresponding section of [1]. Note that to reproduce (exactly) the procedure described in Numerical Recipes [1] one must have npoint=2$\times$M where M is the variable used in the procedure spctrm(), and the decay time must be very large (so that the two successive spectra are equally weighted). If the data being passed is not continuous from one call to the next, set overlap=0.

One frequently wants to do a moving-time average of power spectra, for example to see how the noise spectral properties of an interferometer are changing with time. This is accomplished in avg_spec() by averaging the spectrum with an exponentially-decaying average. Let $A_{t}(f)$ denote the average power spectrum as a function of frequency $f$, at time $t$. Then the exponentially-decaying average ${\langle A(f) \rangle}_t$ at time $t$ is defined by

$${\langle A(f) \rangle}_t = { \int_{-\infty}^t dt' \; A_{t'}... ...u} \over \int_{-\infty}^t dt' \; {\rm e}^{-(t-t')/\tau} },$$ (16.3.321)

where $\tau$ is the characteristic decay time over which an impulse in the power spectrum would decay. In our case, we wish to average the power spectra obtained in the nth pass through the averaging routine. The discrete analog of the previous equation () is

$${\langle A(f) \rangle}_N = { \displaystyle \sum_{n=0}^N A_n... ...\over \displaystyle \sum_{n=0}^N {\rm e}^{-\alpha (N-n) } }.$$ (16.3.322)

Here,

$$\alpha = { {\tt npoint} \over {\tt srate} \times {\tt decaytime} }$$ (16.3.323)

is determined by the averaging time desired. The average defined by () can be easily determined by a recursion relation. We denote the the normalization factor by

$${\cal N}_N = \sum_{n=0}^N {\rm e}^{-\alpha (N-n) }.$$ (16.3.324)

It obeys the (stable) recursion relation ${\cal N}_{N} = 1+ {\rm e}^{-\alpha} {\cal N}_{N-1}$ together with the initial condition ${\cal N}_{-1}=0$. The exponentially-decaying average then satisfies the (stable) recursion relation

$${\langle A(f) \rangle}_{N} = {\rm e}^{-\alpha} {{\cal N}_{N-... ...A_{N}(f) \over {\cal N}_{N} } \quad {\rm for\ } N=0,1,2,\cdots$$ (16.3.325)

(no initial condition is needed). The routine avg_spec() computes the exponentially decaying average by implementing these recursion relations for ${\langle A(f) \rangle}_{N}$ and ${\cal N}_N$.

The units of the output array average[ ] are the square of the units of the input array data[ ] per Hz, i.e.

$${\rm units} \left( {\tt average[\;]} \right) = \left( {\rm units}\left( {\tt data[\;]} \right) \right)^2/\rm Hz.$$ (16.3.326)

The example program calibrate described earlier makes use of the routine avg_spec().

Authors: Bruce Allen, ballen@dirac.phys.uwm.edu and Patrick Brady patrick@tapir.caltech.edu

Comments: See comments for calibrate. Warning: If overlap is turned on, and you pass avg_spec() sets of points that are not continuous, you will introduce discontinous jumps between the data sets, and add lots of peculiar high-frequency garbage to the spectrum.