[R] power in a specific frequency band

[Gabor Grothendieck]

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