Re: [R] dpss.taper for spectral estimation

From: Karim Rahim <karim_at_mast.queensu.ca>
Date: Tue 18 Jul 2006 - 08:05:59 EST


Hi Rouyer,

You can redefine dpss.taper as follows

dpss.taper.2 <- function (n, k, nw = 4, nmax = 2^(ceiling(log(n, 2)))) {

     if (n > nmax)
         stop("length of taper is greater than nmax")
     w <- nw/n
     if (w > 0.5)
         stop("half-bandwidth parameter (w) is greater than 1/2")
     if (k <= 0)
         stop("positive dpss order (k) required")
     v <- matrix(0, nrow = nmax, ncol = (k + 1))

     storage.mode(v) <- "double"
     out <- .Fortran("dpss", nmax = as.integer(nmax), kmax = as.integer(k),
         n = as.integer(n), w = as.double(w), v = v, sig = double(k +
             1), totit = integer(1), sines = double(n), vold = double(n),
         u = double(n), scr1 = double(n), ifault = integer(1),
         PACKAGE = "waveslim")
     return(list(  v=out$v[1:n, 1:k],
                   eigen=1+out$sig[1:k],
                   iter=out$totit,
                   n=n,
                   w=w,
                   ifault=out$ifault) );
}

or you can calculate the eigenvalues from the tapers as done in the tridiagonal dpss calculation at

http://lib.stat.cmu.edu/sapaclisp/multitaper.lisp

Karim



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Tue Jul 18 08:12:08 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Tue 18 Jul 2006 - 10:22:54 EST.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.