empiricalspectrum: Compute the "empirical" spectrum of a time series.

In bspec: Bayesian Spectral Inference

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] * (deltat/N) * abs(DFT(x)(f[j]))^2

where j=0,...,N/2+1, and f[j]=j / (N*deltat) are the Fourier frequencies. deltat 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")

bspec documentation built on April 20, 2022, 5:08 p.m.