### 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 sizeFs <- 128# time vectort  <- seq(0,1,,Fs)# randomize the offsetrt  <- t + rnorm(1)# amplitude of signal 1A1 <- 4# amplitude of signal 2A2 <- 2# annual signal frequencyF1 <- Fs/12# semi annual frequencyF2 <- Fs/6# calculation of yy  <- A1*sin(2*pi*rt*F1) + A2*sin(2*pi*rt*F2)# calculation of powerpower <- (A1^2/2) + (A2^2/2)start_month <- 36# this is the graphmonth <- 1 : Fsplot (month, y, type='l')png("power_calc.png", bg="transparent", width=600, height=600)plot (month, y, type='l')dev.off()# this is the calculationdate()outmat<-NULLfor (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 powerpower# this is what was calculatedoutmat