# setwelch: Set up Matrix of fft for Welch method In RSEIS: Seismic Time Series Analysis Tools

## Description

Prepares a matrix for estimation of power spectrum via Welch's method. Also, is can be used for spectrogram.

## Usage

 ```1 2``` ```setwelch(X, win = min(80, floor(length(X)/10)), inc = min(24, floor(length(X)/30)), coef = 64, wintaper=0.05) ```

## Arguments

 `X` Time series vector `win` window length `inc` increment `coef` coefficient for fft `wintaper` percent taper window taper

## Value

List:

 `values` Matrix of fft's staggered along the trace `windowsize` window length used `increment` increment used `wintaper` percent taper window taper

## Author(s)

originally written by Andreas Weingessel, modified Jonathan M. Lees<[email protected]>

## References

Welch, P.D. (1967) The use of Fast Fourier Transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms IEEE Trans. Audio Electroacoustics 15, 70-73.

stft

## 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64``` ``` dt <- 0.001 t <- seq(0, 6, by=dt) x <- 6*sin(2*pi*50*t) + 10* sin(2*pi*120*t) y <- x + rnorm(length(x), mean=0, sd=10) plot(t,y, type='l') title('sin(2*pi*50*t) + sin(2*pi*120*t)+ rnorm') Y <- fft(y) Pyy <- Y * Conj(Y) N <- length(y) n <- length(Pyy)/2 Syy <- (Mod(Pyy[1:n])^2)/N fn <- 1/(2*dt) f <- (0:(length(Syy)-1))*fn/length(Syy) plot(f, Syy, type='l', log='y' , xlim=c(0, 150)); abline(v=c(50, 120),col='blue', lty=2) plot(f, Syy, type='l', log='y' , xlim=c(0, 150)); abline(v=c(50, 120),col='blue', lty=2) win <- 1024 inc <- min(24, floor(length(y)/30)) coef <- 2048 w <- setwelch(y, win=win, inc=inc, coef=coef, wintaper=0.2) KK <- apply(w\$values, 2, FUN="mean") fw <- seq(from=0, to=0.5, length=coef)/(dt) plot(fw, KK^2, log='', type='l' , xlim=c(0, 150)) ; abline(v=c(50, 120), col='blue', lty=2) Wyy <- (KK^2)/w\$windowsize plot(f, Syy, type='l', log='y' , xlim=c(0, 150)) lines(fw,Wyy , col='red') DBSYY <- 20*log10(Syy/max(Syy)) DBKK <- 20*log10(Wyy/max(Wyy)) plot(f, DBSYY, type='l' , xlim=c(0, 150), ylab="Db", xlab="Hz") lines(fw, DBKK, col='red') title("Compare simple periodogam with Welch's Method") ```

RSEIS documentation built on May 29, 2017, 10:38 p.m.