Estou coletando dados de temperatura de uma geladeira. Os dados parecem uma onda. Gostaria de determinar o período e a frequência da onda (para que eu possa medir se as modificações no refrigerador têm algum efeito).
Estou usando R e acho que preciso usar uma FFT nos dados, mas não tenho certeza para onde ir a partir daí. Eu sou muito novo em análise de sinal e R, então qualquer ajuda seria muito apreciada!
Aqui está a onda que estou produzindo:
Aqui está o meu código R até agora:
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))
Eu publiquei o código R junto com o banco de dados SQLite aqui .
Aqui está um gráfico do gráfico normalizado (com a média removida):
Por enquanto, tudo bem. Aqui está o gráfico de densidade espectral:
Em seguida, aumentamos o zoom no lado esquerdo do gráfico e marcamos o índice mais alto (que é 70) com uma linha verde:
Finalmente, calculamos a frequência da onda. Essa onda é muito lenta, então a convertemos em minutos por ciclo e imprimimos esse valor que é 17,14286.
Aqui estão meus dados em formato delimitado por tabulação, se alguém quiser tentar.
Obrigado pela ajuda! Este problema foi difícil (para mim), mas eu me diverti muito!
Respostas:
Projeto interessante que você está fazendo lá! :-)
De um POV de análise de sinal, essa é realmente uma pergunta simples - e sim, você está certo em utilizar a FFT para esse problema de estimativa de frequência.
Então, muito simplesmente, encontre o máximo de onde está o seu PSD. A abscissa desse máximo corresponderá à sua frequência.
Advertência Emptor, estou dando uma visão geral e suspeito que o resultado da FFT em R seja a frequência normalizada, caso em que você precisaria conhecer sua taxa de amostragem (o que você faz) para convertê-la novamente em Hz. Há muitos outros detalhes importantes que estou deixando de fora, como a resolução da frequência, o tamanho da FFT e o fato de que você provavelmente deseja reduzir o sinal primeiro, mas será bom ver um gráfico primeiro.
EDITAR:
Vamos levar seu sinal em consideração. Depois que eu quero dizer isso, fica assim:
Você pode ver como é simétrico. Se você ignorar a última metade e apenas olhar para a primeira metade e ampliar o zoom, poderá ver o seguinte:
fonte
Para uma forma de onda tão suave e estacionária, a contagem de pontos de amostra entre cruzamentos positivos positivos de algum valor limite médio fornecerá uma estimativa de período. Observe vários períodos de passagem do limite para obter uma estimativa mais média ou detectar qualquer tendência.
fonte
Não há necessidade de fazer nada complicado: basta medir a duração entre os picos da forma de onda. Este é o período. A frequência é apenas 1 dividido pelo período.
Com cerca de 8 ciclos em 2 horas, a frequência é de 4 ciclos por hora, ou cerca de 1 mHz.
fonte