Determinando a frequência e o período de uma onda

11

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:

Minha onda

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):

gráfico normalizado

Por enquanto, tudo bem. Aqui está o gráfico de densidade espectral:

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:

ampliar o gráfico espectral

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!

Aaron Patterson
fonte
Aaron, acho que a melhor coisa aqui é colocar um link para seu arquivo de dados (como texto ou algo assim) em uma caixa de depósito, para que eu possa fazer o download e dar a resposta. Caso contrário, haverá muitas coisas indo e voltando. Não consigo distinguir os números no extremo esquerdo. :-) (Também dê a você a taxa de amostragem - ou seja, com que frequência você faz uma leitura de temperatura).
Spacey
Ah desculpa. Os dados contêm temperatura em graus C, converti em graus F para o gráfico. Porém, são os dados corretos (é a coluna "temp").
Aaron Patterson
O problema com a medição da frequência dessa maneira é que, se você obtiver uma variabilidade considerável de ciclo para ciclo, será muito mais difícil determinar a frequência média - os picos se espalharão juntos - enquanto a simples contagem do tempo entre as excursões permitirá que você calcule bem as coisas (e também calcule std dev, etc). Usar a abordagem FFT seria mais necessário se houvesse muito ruído, mas esse não parece ser o caso aqui.
Daniel R Hicks
+1 para postagem, código, dados, plotagens e um link para o github.
Nibot
@DanielRHicks Nesse caso em particular, não acho que isso importe, mas sim, a FFT fornecerá a média de todas elas, enquanto que, se fizéssemos algo como um cruzamento de zero, mediríamos a duração de cada ciclo (frequência), e então podemos determinar se queremos ter média, mediana, modo etc. Bom ponto!
Spacey

Respostas:

7

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.

real2+imag2

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:

insira a descrição da imagem aqui

x[n]

Nfft=103600=36000.

fs=0.5Hz

x[n]X(f)10log10(|X(f)|2)

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

FsNfft=1.3889e005Hz01.3889e005=0Hz11.3889e005=1.3889e005Hz701.3889e005=9.7222e004Hz

1(9.7222e004)60=17.14

Spacey
fonte
@AaronPatterson Eu editei a publicação, consulte. Além disso, você pode adicionar suas fotos diretamente à sua postagem original. :-). Adicione uma imagem do resultado do PSD que você obtiver.
Spacey
1
Não é exatamente correto se a frequência estiver entre os compartimentos de resultados da FFT.
precisa
@ hotpaw2 Foi por isso que avisei o OP que estou dando uma perspectiva geral e por que preciso ver o enredo. Mesmo assim, editei para adicionar as advertências extras.
Spacey
1
@AaronPatterson Sem problemas, prazer em ajudar. Quanto aos livros, veja Richard Lyons "Understanding DSP" - este é um livro rápido para começar.
Spacey
1
1.3x105
4

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.

hotpaw2
fonte
3

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.

nibot
fonte
3
Como posso fazer isso programaticamente?
Aaron Patterson