Unfortunately the documentation seems somewhat sparse
but you can figure out what it is doing by looking at the code
of spec.pgram or just pumping a sine wave through it and
comparing that to a manually generated periodogram, which
we do below:
> tt <- 1:32 # time
> x <- sin(2*pi*tt/8) # sine save with period 8, i.e. freq = .125
> x <- x - mean(x)
>
> # Calculate the periodogram manually using the fft
> round(Mod(fft(x))^2/32, 3)
[1] 0 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 0
0>
> # Now recalculate using spectrum
> sp <- spectrum(x, taper = 0, detrend = FALSE)
>
> # Note the power concentrated at sp$freq[4]=.125 in sp$spec[4].
> # Note that the fft version returns the mean in element 1 so element 2
> # of the fft version corresponds to element 1 of the spectrum version.
> # Also in comparing we see that the spectrum version omits the symmetric
> # completion of the first half.
> round(sp$spec,2)
[1] 0 0 0 8 0 0 0 0 0 0 0 0 0 0 0 0> sp$freq
[1] 0.03125 0.06250 0.09375 0.12500 0.15625 0.18750 0.21875 0.25000 0.28125
[10] 0.31250 0.34375 0.37500 0.40625 0.43750 0.46875
0.50000> 1/sp$freq
[1] 32.000000 16.000000 10.666667 8.000000 6.400000 5.333333 4.571429
[8] 4.000000 3.555556 3.200000 2.909091 2.666667 2.461538 2.285714
[15] 2.133333 2.000000>
> # sum of squares vs. sum of periodogram values.
> # Note that power defined as sum of squares of demeaned series
> # equals twice sum of periodogram values returned by spectrum
> # (since spec.pgram does not return the half that is the same
> # by symmetry)
> sum(x*x)
[1] 16> sum(sp$spec)
[1] 8
Date: Fri, 15 Oct 2004 19:55:40 -0400
From: Giuseppe Pagnoni <gpagnon at emory.edu>
To: <r-help at stat.math.ethz.ch>
Subject: [R] power in a specific frequency band
Dear R users
I have a really simple question (hoping for a really simple answer :-):
Having estimated the spectral density of a time series "x" (heart rate
data) with:
x.pgram <- spectrum(x,method="pgram")
I would like to compute the power in a specific energy band.
Assuming that frequency(x)=4 (Hz), and that I am interested in the band
between f1 and f2, is the power in the band simply the following?
sum(x.pgram$spec[(x.pgram$freq > f1/frequency(x)) & (x.pgram$freq <=
f2/frequency(x))])
If it is so, are the returned units the same units as the original time
series, but squared (if x is bpm, then the power is in bpm^2)?
I own a copy of Venables and Ripley (MASS 2003), but I was not able to
extract this information from the time series chapter....
thanks for any help
(please cc to my e-mail, if possible)
giuseppe
Giuseppe Pagnoni
Dept. Psychiatry and Behavioral Sciences
Emory University
1639 Pierce Drive, Suite 4000
Atlanta, GA, 30322
tel: 404.712.8431
fax: 404.727.3233