El uso de una ventana móvil.
Por ejemplo, considere estos datos sintéticos (creado en R):
time <- 1:(365*24)/24
amplitude <- time * ((1+length(time))/24-time) / length(time)
set.seed(17)
x <- (x0 <- amplitude * cos(2*pi * (time + 6 * (time/365)^2))) + rnorm(length(time), sd=5)
![Time series]()
Los puntos oscuros son los que subyacen silencioso de la señal; la línea gris es la señal con el agregado de ruido Gaussiano (croma en 7 intervalos de hora para hacerlo más claro). La amplitud de la señal obviamente, pero un poco menos obviamente, la frecuencia aumenta con el tiempo, demasiado.
Ejecución de una ventana a través de esta serie, la informática, el periodograma para cada posición, y extraer el valor de la frecuencia dominante encuentra en la primera ventana, produce una serie de facultades que corresponden a los centros de windows. Si no hay ningún cambio en el poder a lo largo del tiempo, la trama de esta serie debería ser casi plana. Las desviaciones pueden ser interpretados como cambios temporales en el poder. Es útil para suavizar esta trama; un lowess suave (con una estrecha ventana de 1/4 de la anchura total) se utiliza aquí.
![Results]()
Como base para la comparación, el mismo procedimiento se aplicó a la base de la señal; se traza en gris. El alisado de la ventana de energía de la función se aproxima a la ventana de la función de potencia para el subyacente de la señal (y está muy cerca en el mayor de los poderes).
(Los poderes que se muestran aquí son calculadas usando logaritmos naturales, no como los decibelios.)
Uno podría hacer un complot para cualquier frecuencia deseada, además de esta "dominante" o la frecuencia esperada.
Apéndice
Aquí es el R código que genera las figuras.
#
# Return the square of the periodogram (normalized by the length of x)
# evaluated at index n.
#
power <- function(x, n) {
x.hat <- spec.pgram(x, plot=FALSE)
(x.hat$spec[n] / length(x))^2 #$(TeX bug workaround)
}
#
# Window `power(*,n)` across array `x` using a weighted symmetric
# window extending to `width` on either side (and therefore of length
# 2*width+1), subsampled every `skip` elements.
#
power.window <- function(x, n, width=1, weight=1, skip=1) {
i <- seq.int(from=width+1, to=length(x)-width, by=skip)
sapply(i, function(i) power(x[(i-width):(i+width)] * weight, n))
}
#
# Simulate and plot an interesting time series.
#
time <- 1:(365*24)/24
amplitude <- time * ((1+length(time))/24-time) / length(time)
set.seed(17)
x0 <- amplitude * cos(2*pi * (time + 6 * (time/365)^2))
x <- x0 + rnorm(length(time), sd=5)
i <- seq.int(from=1, to=length(time), by=7)
plot(time[i], x[i], type="l", xlab="Days", ylab="x", main="Data", col="Gray")
points(time, x0, pch=19, cex=0.2)
#
# Find the frequency where power is maximized within the initial window.
#
width <- 30*24
x.hat <- spec.pgram(x[1:(2*width+1)], plot=FALSE)
i.max <- which.max(x.hat$spec) #$ TeX bug workaround
#
# Compute and plot the power moving window for the underlying
# series (`x0`) and its noisy version (`x`). For plotting, subsample
# the series.
#
skip <- 7
x.power <- power.window(x, n=i.max, width=width, weight=1, skip=skip)
x0.power <- power.window(x0, n=i.max, width=width, weight=1, skip=skip)
x.power.smooth <- lowess(x.power, f=1/4)
i <- seq(from=(width+1)/24, by=skip/24, length.out=length(x.power))
main <- sprintf("Moving power at 1/%4.1f hours, %d hour window",
1/x.hat$freq[i.max], 2*width+1) #$ TeX bug
plot(i, log(x.power), type="l", lwd=2, xlab="Time (days)", main=main)
lines(i, log(x.power.smooth$y), lty=2, lwd=3, col="Blue")
lines(i, log(x0.power), lwd=2, col="Gray")
legend(30, -4, "Underlying", bty="n", box.col="White", box.lwd=0, lwd=2, col="Gray")
legend(30, -5, "Windowed", bty="n", box.col="White", box.lwd=0, lwd=2)
legend(30, -6, "Smooth", bty="n", box.col="White", box.lwd=0, lwd=3, col="Blue", lty=2)
