power_demo

> #
> # Finding a signal's power using different methods
> #
> # This was mostly cobbled together using the examples from the documentation
> # for the pwelch function from the R package "oce" (oceanography)
> #
> # http://finzi.psych.upenn.edu/R/library/oce/html/pwelch.html
> #
> #   (Thank-you Dan Kelley!... I owe you a Tim Horton's)
> # 
> # I also used an example from the pwelch documentation from MATLAB, but
> # I can't seem to find it anymore.
> #
> library(oce)
> # sample size
> Fs <- 128
> # time vector
> t  <- seq(0,1,,Fs)
> # randomize the offset
> rt  <- t + rnorm(1)
> # amplitude of signal 1
> A1 <- 4
> # amplitude of signal 2
> A2 <- 2
> # annual signal frequency
> F1 <- Fs/12
> # semi annual frequency
> F2 <- Fs/6
> # calculation of y
> y  <- A1*sin(2*pi*rt*F1) + A2*sin(2*pi*rt*F2)
> # calculation of power
> power <- (A1^2/2) + (A2^2/2)
> start_month <- 36
> # this is the graph
> month <- 1 : Fs
> plot (month, y, type='l')



> png("power_calc.png", bg="transparent", width=600, height=600)
> plot (month, y, type='l')
> dev.off()
windows 
      2 
> # this is the calculation
> date()
[1] "Sat Jul 30 04:04:02 2011"
> outmat<-NULL
> for (i in start_month:128){
+   yi <- y[1:i]
+   yts <- ts(yi, frequency=Fs)
+   
+   w1 <- pwelch(yts)
+   w2 <- pwelch(yts, nfft=48)
+   
+   pvar <- var(yi)
+   pabsq <- sum(abs(yi)^2)/length(yi)
+   pwelch <- sum(w1$spec) * diff(w1$freq[1:2])
+   pwelch_nfft <- sum(w2$spec) * diff(w2$freq[1:2])
+   
+   month <- i
+   
+   out<-cbind(month, pvar,pabsq,pwelch,pwelch_nfft)
+   
+   outmat<-rbind(outmat,out)
+ }
> date()
[1] "Sat Jul 30 04:04:15 2011"
> # this is the power
> power
[1] 10
> # this is what was calculated
> outmat
      month      pvar     pabsq   pwelch pwelch_nfft
 [1,]    36 10.379178 10.092213 5.274233    9.992992
 [2,]    37 10.810329 10.549180 5.361375   10.010590
 [3,]    38 10.899800 10.689332 5.361375   10.013920
 [4,]    39 10.664935 10.489352 5.361375   10.007797
 [5,]    40 10.391572 10.228687 2.511616    9.999187
 [6,]    41 10.134147  9.979207 2.449863    9.991234
 [7,]    42  9.895617  9.743733 2.449863    9.985259
 [8,]    43  9.764578  9.595115 2.449863    9.982640
 [9,]    44  9.968130  9.761493 2.449863    9.986311
[10,]    45 10.373281 10.143276 2.418762    9.993470
[11,]    46 10.416288 10.192815 2.418762    9.995721
[12,]    47 10.192957  9.978234 2.418762    9.995178
[13,]    48 10.312052 10.098610 3.440111   10.000055
[14,]    49 10.637727 10.440892 3.440111   10.006895
[15,]    50 10.685453 10.518020 3.664137   10.007745
[16,]    51 10.504522 10.356378 3.664137   10.004798
[17,]    52 10.298606 10.157888 3.664137   10.001428
[18,]    53 10.101643  9.966230 3.664137    9.998321
[19,]    54  9.918106  9.784380 3.441609    9.995536
[20,]    55  9.827327  9.681911 3.441609    9.994046
[21,]    56 10.007271  9.839047 6.225073    9.996495
[22,]    57 10.314496 10.133640 5.718061   10.000183
[23,]    58 10.314328 10.138593 5.718061    9.999531
[24,]    59 10.146472  9.975574 5.718061    9.997577
[25,]    60 10.274173 10.104371 5.718061   10.000800
[26,]    61 10.530787 10.372943 5.718061   10.005600
[27,]    62 10.549983 10.411239 5.771936   10.005920
[28,]    63 10.400764 10.273884 5.771936   10.003698
[29,]    64 10.235731 10.113637 6.239504   10.001359
[30,]    65 10.076398  9.958043 6.239504    9.999183
[31,]    66  9.928185  9.810638 6.239504    9.997140
[32,]    67  9.865912  9.739701 6.239504    9.996162
[33,]    68 10.032194  9.890614 5.775204    9.998365
[34,]    69 10.273728 10.124838 5.775204   10.001086
[35,]    70 10.246197 10.101389 5.775204    9.999936
[36,]    71 10.118455  9.976488 5.775204    9.998265
[37,]    72 10.250344 10.109449 4.314942   10.001118
[38,]    73 10.456984 10.325295 4.314942   10.004812
[39,]    74 10.456144 10.337758 4.314942   10.004879
[40,]    75 10.328259 10.217689 4.718504   10.002987
[41,]    76 10.190617 10.083359 4.718504   10.001082
[42,]    77 10.056918  9.952407 4.718504    9.999271
[43,]    78  9.933513  9.829246 4.718504    9.997564
[44,]    79  9.892925  9.781812 4.718504    9.996910
[45,]    80 10.049708  9.927618 5.688715    9.999019
[46,]    81 10.242745 10.116302 5.688715   10.001184
[47,]    82 10.197242 10.074091 5.350431    9.999907
[48,]    83 10.101004  9.979571 5.350431    9.998660
[49,]    84 10.234142 10.113810 5.350431   10.001343
[50,]    85 10.402198 10.289262 5.350431   10.004358
[51,]    86 10.387010 10.283780 5.350431   10.004147
[52,]    87 10.274771 10.176943 5.350431   10.002487
[53,]    88 10.156747 10.061337 7.408198   10.000821
[54,]    89 10.041620  9.948292 7.408198    9.999229
[55,]    90  9.936820  9.843362 7.408198    9.997747
[56,]    91  9.913612  9.814519 7.408198    9.997346
[57,]    92 10.062729  9.955494 7.408198    9.999384
[58,]    93 10.217620 10.107802 7.408198   10.001132
[59,]    94 10.160386 10.053238 7.408198    9.999820
[60,]    95 10.090175  9.984077 7.089592    9.998968
[61,]    96 10.222381 10.117428 7.646170   10.000870
[62,]    97 10.359330 10.260491 7.646170   10.000870
[63,]    98 10.333786 10.242263 7.646170   10.000169
[64,]    99 10.233701 10.146046 7.646170   10.000169
[65,]   100 10.130410 10.044599 7.646170   10.000169
[66,]   101 10.029361  9.945160 7.646170   10.000169
[67,]   102  9.939252  9.854666 7.646170   10.000169
[68,]   103  9.930504  9.841154 7.752149   10.000169
[69,]   104 10.072697  9.977151 8.143510   10.000169
[70,]   105 10.196270 10.099256 8.143510   10.000169
[71,]   106 10.131792 10.036950 8.143510   10.000169
[72,]   107 10.083752  9.989549 8.143510   10.000169
[73,]   108 10.213308 10.120291 8.143510   10.000169
[74,]   109 10.324435 10.236575 8.123304   10.000169
[75,]   110 10.291447 10.209234 8.123304   10.000169
[76,]   111 10.201185 10.121818 8.123304   10.000169
[77,]   112 10.109359 10.031448 8.599350   10.000169
[78,]   113 10.019355  9.942706 8.599350   10.000169
[79,]   114  9.941355  9.864141 8.599350   10.000169
[80,]   115  9.944949  9.863634 8.599350   10.000169
[81,]   116 10.080394  9.994284 8.599350   10.000169
[82,]   117 10.177509 10.090655 8.599350   10.000169
[83,]   118 10.109202 10.024113 8.279687   10.000169
[84,]   119 10.080377  9.995673 8.279687   10.000169
[85,]   120 10.205880 10.122398 8.816538   10.000169
[86,]   121 10.295155 10.216084 8.816538   10.000169
[87,]   122 10.256912 10.182274 8.816538   10.000169
[88,]   123 10.174807 10.102313 8.962555    9.999720
[89,]   124 10.092153 10.020843 8.962555    9.999720
[90,]   125 10.011062  9.940744 8.962555    9.999720
[91,]   126  9.943410  9.872400 8.962555    9.999720
[92,]   127  9.957713  9.883128 8.962555    9.999720
[93,]   128 10.086278 10.007946 9.275188    9.999720

XXXXXXXXXXXXXXXXXXXX  SOURCE CODE XXXXXXXXXXXXXXXXX

#
# Finding a signal's power using different methods
#
# This was mostly cobbled together using the examples from the documentation
# for the pwelch function from the R package "oce" (oceanography)
#
# http://finzi.psych.upenn.edu/R/library/oce/html/pwelch.html
#
#   (Thank-you Dan Kelley!... I owe you a Tim Horton's)
# I also used an example from the pwelch documentation from MATLAB, but
# I can't seem to find it anymore.
#

library(oce)

# sample size
Fs <- 128

# time vector
t  <- seq(0,1,,Fs)

# randomize the offset
rt  <- t + rnorm(1)

# amplitude of signal 1
A1 <- 4

# amplitude of signal 2
A2 <- 2

# annual signal frequency
F1 <- Fs/12

# semi annual frequency
F2 <- Fs/6

# calculation of y
y  <- A1*sin(2*pi*rt*F1) + A2*sin(2*pi*rt*F2)

# calculation of power
power <- (A1^2/2) + (A2^2/2)

start_month <- 36

# this is the graph

month <- 1 : Fs

plot (month, y, type='l')

png("power_calc.png", bg="transparent", width=600, height=600)

plot (month, y, type='l')

dev.off()

# this is the calculation

date()
outmat<-NULL

for (i in start_month:128){

  yi <- y[1:i]
  yts <- ts(yi, frequency=Fs)
  
  w1 <- pwelch(yts)
  w2 <- pwelch(yts, nfft=48)
  
  pvar <- var(yi)
  pabsq <- sum(abs(yi)^2)/length(yi)
  pwelch <- sum(w1$spec) * diff(w1$freq[1:2])
  pwelch_nfft <- sum(w2$spec) * diff(w2$freq[1:2])
  
  month <- i
  
  out<-cbind(month, pvar,pabsq,pwelch,pwelch_nfft)
  
  outmat<-rbind(outmat,out)
}
date()

# this is the power
power

# this is what was calculated
outmat

Comments