# [R] power in a specific frequency band

Gabor Grothendieck ggrothendieck at myway.com
Sat Oct 16 19:16:48 CEST 2004

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

```