btpsd: Blackman-Tukey spectral density estimate

Description Arguments Value Note Author(s) References See Also Examples

View source: R/BTMEM.R

Description

Computes the Blackman-Tukey spectral density estimate

Arguments

y

time series.

type

Character. Default is "Tukey". Other window choices are "Triangular","Hanning",and "Hamming".

win

window length.

taper

amount to taper the time series. Default is 0.5.

Value

est

Blackman-Tukey estimates of the power spectral density.

Note

~~further notes~~

Author(s)

Patrick Crutcher

References

D.B. Percival and A. Walden (1993), Spectral Analysis for Physical Applications: Multitaper and Conventional Univariate Techniques, Cambridge University Press, Cambridge, UK.

See Also

~~objects to See Also as help, ~~~

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function(y, type="Tukey", win=N , taper=0.5) {
  
  y <- y - mean(y)
  T <- length(y)
  if ( win >= T ) stop("The length of the window is longer than the data length.")

  if (is.null(win)) N <- ceiling(2*sqrt(T)) else N <- win
  
  if (type =="Tukey") { w <- tukey(win,taper) }
  if (type =="Hanning") { w <- hanning.window(win) }
  if (type =="Hamming") { w <- hamming.window(win) }
  if (type =="Triangular") { w <- tri.window(win) }


	r <- acf(y, lag.max = N-1,plot=FALSE)$acf;

	rw <- r*w[-N]
	
	est <- Mod(fft(rw))^2/(2*pi*length(y))
	return(est)
}

Example output

function (y, type = "Tukey", win = N, taper = 0.5) 
{
    y <- y - mean(y)
    T <- length(y)
    if (win >= T) 
        stop("The length of the window is longer than the data length.")
    if (is.null(win)) 
        N <- ceiling(2 * sqrt(T))
    else N <- win
    if (type == "Tukey") {
        w <- tukey(win, taper)
    }
    if (type == "Hanning") {
        w <- hanning.window(win)
    }
    if (type == "Hamming") {
        w <- hamming.window(win)
    }
    if (type == "Triangular") {
        w <- tri.window(win)
    }
    r <- acf(y, lag.max = N - 1, plot = FALSE)$acf
    rw <- r * w[-N]
    est <- Mod(fft(rw))^2/(2 * pi * length(y))
    return(est)
}

ssa documentation built on May 2, 2019, 5:54 p.m.

Related to btpsd in ssa...