Determinación de frecuencia y período de una ola

11

Estoy recopilando datos de temperatura de un refrigerador. Los datos se ven como una ola. Me gustaría determinar el período y la frecuencia de la ola (para poder medir si las modificaciones en el refrigerador tienen algún efecto).

Estoy usando R, y creo que necesito usar una FFT en los datos, pero no estoy seguro de dónde ir desde allí. Soy muy nuevo en R y análisis de señales, por lo que cualquier ayuda sería muy apreciada.

Aquí está la ola que estoy produciendo:

Mi ola

Aquí está mi código R hasta ahora:

require(graphics)
library(DBI)
library(RSQLite)

drv <- dbDriver("SQLite")
conn <- dbConnect(drv, dbname = "s.sqlite3")

query <- function(con, query) {
  rs <- dbSendQuery(con, query)
  data <- fetch(rs, n = -1)
  dbClearResult(rs)
  data
}

box <- query(conn, "
SELECT id,
       humidity / 10.0 as humidity,
       temp / 10.0 as temp,
       ambient_temp / 10.0 as ambient_temp,
       ambient_humidity / 10.0 as ambient_humidity,
       created_at
FROM measurements ORDER BY id DESC LIMIT 3600
")

box$x <- as.POSIXct(box$created_at, tz = "UTC")

box$x_n <- box$temp - mean(box$temp)
png(filename = "normalized.png", height = 750, width = 1000, bg = "white")
plot(box$x, box$x_n, type="l")

# Pad the de-meaned signal so the length is 10 * 3600
N_fft  <- 3600 * 10
padded <- c(box$x_n, seq(0, 0, length= (N_fft - length(box$x_n))))
X_f    <- fft(padded)
PSD    <- 10 * log10(abs(X_f) ** 2)

png(filename = "PSD.png", height = 750, width = 1000, bg = "white")
plot(PSD, type="line")

zoom <- PSD[1:300]

png(filename = "zoom.png", height = 750, width = 1000, bg = "white")
plot(zoom, type="l")

# Find the index with the highest point on the left half
index <- which(PSD == max(PSD[1:length(PSD) / 2]))

# Mark it in green on the zoomed in graph
abline(v = index, col="green")

f_s     <- 0.5 # sample rate in Hz
wave_hz <- index * (f_s / N_fft)
print(1 / (wave_hz * 60))

He publicado el código R junto con la base de datos SQLite aquí .

Aquí hay una gráfica del gráfico normalizado (con la media eliminada):

gráfico normalizado

Hasta aquí todo bien. Aquí está la gráfica de densidad espectral:

densidad espectral

Luego, ampliamos el lado izquierdo del gráfico y marcamos el índice más alto (que es 70) con una línea verde:

acercar el diagrama espectral

Finalmente calculamos la frecuencia de la onda. Esta onda es muy lenta, por lo que la convertimos a minutos por ciclo e imprimimos ese valor que es 17.14286.

Aquí están mis datos en formato delimitado por tabuladores si alguien más quiere probar.

¡Gracias por la ayuda! ¡Este problema fue difícil (para mí) pero lo pasé muy bien!

Aaron Patterson
fuente
Aaron, creo que lo mejor aquí es que pongas un enlace a tu archivo de datos (como texto o algo) en un buzón, para que pueda descargarlo y darte la respuesta. De lo contrario, habrá muchos cambios de ida y vuelta. No puedo distinguir los números en el extremo izquierdo. :-) (También le da una frecuencia de muestreo, es decir, con qué frecuencia está tomando una lectura de temperatura).
Spacey
Oh, lo siento. Los datos contienen temperatura en grados C, convertí a grados F para el gráfico. Sin embargo, son los datos correctos (es la columna "temp").
Aaron Patterson
El problema con la medición de la frecuencia de esa manera es que si obtiene una variabilidad considerable de un ciclo a otro, será mucho más difícil determinar la frecuencia promedio (los picos se unirán), mientras que simplemente contar el tiempo entre excursiones le permitirá promediar las cosas muy bien (y también calcular std dev, etc.) Usar el enfoque FFT sería más necesario si hubiera mucho ruido, pero ese no parece ser el caso aquí.
Daniel R Hicks
+1 para publicación, código, datos, tramas y un enlace a github.
nibot
@DanielRHicks En este caso particular, no creo que importe, pero sí, la FFT le dará el promedio de todos ellos, mientras que si hiciéramos algo como un cruce por cero, mediríamos la duración (frecuencia) de cada ciclo, y luego podemos determinar si queremos tomar la media, la mediana, la moda, etc. ¡Buen punto!
Spacey

Respuestas:

7

¡Interesante proyecto que tienes allí! :-)

Desde un punto de vista de análisis de señal, esta es realmente una pregunta simple, y sí, tiene razón en utilizar la FFT para este problema de estimación de frecuencia.

No estoy familiarizado con R, pero lo que esencialmente quieres hacer es tomar la FFT de tu señal de temperatura. Como su señal es real, obtendrá un resultado complejo como la salida de su FFT. Ahora puede tomar la magnitud absoluta al cuadrado de su resultado FFT complejo (ya que no nos importa la fase). Es decir, . Esta es tu densidad espectral de potencia.real2+imag2

Luego, de manera muy simple, encuentre el máximo de dónde se encuentra su PSD. La abscisa de ese máximo corresponderá a su frecuencia.

Advertencia: Emptor, le estoy dando una perspectiva general, y sospecho que el resultado de la FFT en R será la frecuencia normalizada, en cuyo caso tendría que conocer su frecuencia de muestreo (lo que hace), para poder volver a convertirla. en Hz. Hay muchos otros detalles importantes que estoy omitiendo, como su resolución de frecuencia, el tamaño de FFT y el hecho de que probablemente desee desviar su señal primero, pero será bueno ver primero una trama.

EDITAR:

Permítanos tener en cuenta su señal. Después de decirlo en serio, se ve así:

ingrese la descripción de la imagen aquí

Desea reducir el significado porque la polarización de CC es una frecuencia de 0 Hz técnicamente, y no desea que ahogue otras frecuencias de interés. Llamemos a esta señal desmedida .x[n]

Ahora, cuando toma la FFT, debe prestar atención a algunos detalles. Utilizo una longitud FFT de 10 veces la longitud de la señal, por lo que podemos decir que Cualquier tamaño de FFT establecido mayor que el tamaño de la señal original tiene el efecto de interpolar su resultado FFT en el dominio de la frecuencia. Si bien esto no agrega ninguna información para usted desde un punto de vista de la teoría de la información, lo ayuda a discernir mejor dónde se encuentra exactamente su pico real, especialmente en los casos en que se encuentra entre dos intervalos de frecuencia. Para obtener una mejor resolución de frecuencia, desea obtener más datos. (es decir, deje que el sensor funcione durante un período de tiempo más largo).Nfft=103600=36000.

Avanzando, desde el archivo de texto, veo que su sensor está tomando lecturas de temperatura cada 2 segundos. Esto significa que su frecuencia de muestreo , ofs=0.5Hz

Por lo tanto, tomemos una FFT de , y llamemos a ese resultado . Este resultado será complejo y, de hecho, es simétrico conjugado, pero eso no es importante aquí. (Solo significa para sus propósitos, la mitad es redundante). Ahora, si tomamos , esta es la densidad espectral de potencia (en escala logarítmica). El resultado se ve así:x[n]X(f)10log10(|X(f)|2)

ingrese la descripción de la imagen aquí

Puedes ver cómo es simétrico. Si ignoras la última mitad, y solo miras la primera mitad y te acercas, puedes ver esto:

ingrese la descripción de la imagen aquí

¡Entonces tienes un pico en el índice 70! Pero, ¿qué es el índice 70 en términos de frecuencia real? Aquí es donde quieres tu resolución de frecuencia. Para calcular eso, simplemente tomamos . Esto simplemente significa que cada índice representa una frecuencia que es un múltiplo de este número. El índice 0 es . El índice 1 es . El índice 70 es .0*1,3889e-005=0Hz1*1,3889e-005=1,3889e-005Hz70*1,3889e-005=9,7222e-004HzFsNfft=1.3889e005Hz01.3889e005=0Hz11.3889e005=1.3889e005Hz701.3889e005=9.7222e004Hz

Casi terminamos. Ahora que tiene esta cifra, podemos convertirla en algo más apetecible. Esta cifra solo dice que su señal completa 9.7222e-004 ciclos por segundo. Entonces podemos preguntar, ¿cuántos minutos se necesitan para calcular un ciclo? Por lo tanto, minutos por ciclo.1(9.7222e004)60=17.14

Spacey
fuente
@AaronPatterson He editado la publicación, por favor vea. Además, puede agregar sus imágenes directamente a su publicación original. :-). Agregue una imagen del resultado de PSD que obtenga.
Spacey
1
No es exactamente correcto si la frecuencia resulta estar entre los contenedores de resultados FFT.
hotpaw2
@ hotpaw2 Es por eso que le advertí a OP que estoy dando una perspectiva general, y por qué necesito ver la trama. De todos modos, he editado para agregar las advertencias adicionales.
Spacey
1
@AaronPatterson No hay problema, me alegra poder ayudar. En cuanto a los libros, mira Richard Lyons "Understanding DSP", que es un libro fascinante para comenzar.
Spacey
1
1.3x105
4

Para una forma de onda tan suave y estacionaria, contar puntos de muestra entre cruces positivos en curso de algún valor umbral promedio le dará una estimación del período. Mire varios períodos de cruce de umbral para obtener una estimación más promedio o detectar cualquier tendencia.

hotpaw2
fuente
3

No hay necesidad de hacer nada complicado: solo mide la duración entre los picos de la forma de onda. Este es el periodo. La frecuencia es solo 1 dividida por el período.

Con aproximadamente 8 ciclos durante 2 horas, la frecuencia es de 4 ciclos por hora, o aproximadamente 1 mHz.

nibot
fuente
3
¿Cómo puedo hacer esto mediante programación?
Aaron Patterson