empiricalSpectrum {bspec}

Compute the "empirical" spectrum of a time series.

Description

Computes the "empirical power" of a time series via a discrete Fourier transform.

Usage

empiricalSpectrum(x, two.sided=FALSE)

Arguments

x

a time series (ts object).

two.sided

a logical flag indicating whether the output is supposed to correspond to the one-sided (default) or two-sided spectrum.

Details

Performs a Fourier transform, and then derives (based on the additional information on sampling rate etc. provided via the time series' attributes) the spectral power as a function of frequency. The result is simpler (in a way) than the spectrum() function's output, see also the example below. What is returned is the real-valued frequency series

\kappa_j\frac{\Delta_t}{N}\bigl|\tilde{x}(f_j)\bigr|^2

where j=0,...,N/2+1, and f_j=\frac{j}{N \Delta_t} are the Fourier frequencies. \Delta_t is the time series' sampling interval and N is its length. \kappa_j is =1 for zero and Nyquist frequencies, and =2 otherwise, and denotes the number of (by definition) non-zero Fourier coefficients. In case two.sided=TRUE, the \kappa_j prefactor is omitted.

For actual spectral estimation purposes, the use of a windowing function (see e.g. the tukeywindow() function) is highly recommended.

Value

A list containing the following elements:

frequency

the Fourier frequencies.

power

the spectral power.

kappa

the number of (by definition) non-zero imaginary components of the Fourier series.

two.sided

a logical flag indicating one- or two-sidedness.

Author(s)

Christian Roever, [e-mail]

References

Roever, C., Meyer, R., Christensen, N. Modelling coloured residual noise in gravitational-wave signal processing. Classical and Quantum Gravity, 28(1):015010, 2011. doi: 10.1088/0264-9381/28/1/015010. See also arXiv preprint 0804.3853.

See Also

fft, spectrum, tukeywindow, welchPSD

Examples

# load example data:
data(lh)

# compute spectrum:
spec1 <- empiricalSpectrum(lh)
plot(spec1$frequency, spec1$power, log="y", type="b")

# plot "spectrum()" function's result in comparison:
spec2 <- spectrum(lh, plot=FALSE)
lines(spec2$freq, spec2$spec, col="red")

# make both spectra match:
spec3 <- empiricalSpectrum(lh, two.sided=TRUE)
spec4 <- spectrum(lh, plot=FALSE, taper=0, fast=FALSE, detrend=FALSE)
plot(spec3$frequency, spec3$power, log="y", type="b")
lines(spec4$freq, spec4$spec, col="green")