Mikhail Titov
2009-Jun-19 18:41 UTC
[R] typo in Lomb-Scargle periodogram implementation in spec.ls() from cts package?
Hello! I tried to contact author of the package, but I got no reply. That is why I write it here. This might be useful for those who were using cts for spectral analysis of non-uniformly spaced data. In file spec.ls.R from cts_1.0-1.tar.gz lines 59-60 are written as pgram[k, i, j] <- 0.5 * ((sum(x[1:length(ti)]* cos(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((cos(2 * pi * freq.temp[k] * (ti - tao)))^2) + (sum(x[1:length(ti)] * sin(2 * pi * freq.temp[k] * (ti - tao))))^2 ===> ) <=== /sum((sin(2 * pi * freq.temp[k] * (ti - tao)))^2) Is there a misplaced bracket (shown like ===> ) <===)? Should it be like the following? pgram[k, i, j] <- 0.5 * ((sum(x[1:length(ti)]* cos(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((cos(2 * pi * freq.temp[k] * (ti - tao)))^2) + (sum(x[1:length(ti)] * sin(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((sin(2 * pi * freq.temp[k] * (ti - tao)))^2) ===> ) <== Here is quick reference http://en.wikipedia.org/wiki/Least-squares_spectral_analysis#The_Lomb.E2.80.93Scargle_periodogram . One half coefficient was not applied to entire expression. Also I find weird next lines (61-62) pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, i, j]) First of all, such things should not be in the for loop. Second, I don't quite understand the meaning of it. P.S. Should I use tapering of my data? If I just try to fit sine and cosine, I may not use it, however for FFT windowing is a must. What about Lomb-Scargle? Mikhail
Uwe Ligges
2009-Jun-20 15:54 UTC
[R] typo in Lomb-Scargle periodogram implementation in spec.ls() from cts package?
Please report problems in the code to the package maintainer (CCing). Best, Uwe Ligges Mikhail Titov wrote:> Hello! > > I tried to contact author of the package, but I got no reply. That is why I write it here. This might be useful for those who were using cts for spectral analysis of non-uniformly spaced data. > > In file spec.ls.R from cts_1.0-1.tar.gz lines 59-60 are written as > > pgram[k, i, j] <- 0.5 * ((sum(x[1:length(ti)]* cos(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((cos(2 * > pi * freq.temp[k] * (ti - tao)))^2) + (sum(x[1:length(ti)] * sin(2 * pi * freq.temp[k] * (ti - tao))))^2 ===> ) <=== /sum((sin(2 * pi * freq.temp[k] * (ti - tao)))^2) > > Is there a misplaced bracket (shown like ===> ) <===)? Should it be like the following? > > pgram[k, i, j] <- 0.5 * ((sum(x[1:length(ti)]* cos(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((cos(2 * > pi * freq.temp[k] * (ti - tao)))^2) + (sum(x[1:length(ti)] * sin(2 * pi * freq.temp[k] * (ti - tao))))^2/sum((sin(2 * pi * freq.temp[k] * (ti - tao)))^2) ===> ) <==> > > Here is quick reference http://en.wikipedia.org/wiki/Least-squares_spectral_analysis#The_Lomb.E2.80.93Scargle_periodogram . One half coefficient was not applied to entire expression. > > Also I find weird next lines (61-62) > > pgram[1, i, j] <- 0.5 * (pgram[2, i, j] + pgram[N, i, j]) > > First of all, such things should not be in the for loop. Second, I don't quite understand the meaning of it. > > P.S. Should I use tapering of my data? If I just try to fit sine and cosine, I may not use it, however for FFT windowing is a must. What about Lomb-Scargle? > > Mikhail > > ______________________________________________ > R-help at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code.