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

(16.3.320) `average[N]`. The units of`average[ ]`are . Note: the elements of`average[ ]`must not be changed in between successive calls to`avg_spec()`.`npoint:`Input. The number of points 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 of the input data, in Hz.`decaytime:`Input. The characteristic (positive) decay time 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 .`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 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`=2M 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 `o`verlap=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
denote the average power spectrum as a function of frequency
, at time . Then the exponentially-decaying average
at time is defined by

Here,

(16.3.323) |

(16.3.324) |

(16.3.325) |

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

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.