Zbieram dane o temperaturze z lodówki. Dane wyglądają jak fala. Chciałbym określić okres i częstotliwość fali (aby móc zmierzyć, czy modyfikacje lodówki mają jakikolwiek wpływ).
Używam R i myślę, że muszę użyć FFT na danych, ale nie jestem pewien, dokąd się udać. Jestem bardzo nowy w R i analizie sygnałów, więc każda pomoc byłaby bardzo mile widziana!
Oto fala, którą produkuję:
Oto mój kod R do tej pory:
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))
Mam pisał kod R wraz z bazą danych SQLite tutaj .
Oto wykres znormalizowanego wykresu (z usuniętą średnią):
Na razie w porządku. Oto wykres gęstości widmowej:
Następnie powiększamy po lewej stronie wykresu i zaznaczamy najwyższy wskaźnik (czyli 70) zieloną linią:
Na koniec obliczamy częstotliwość fali. Fala jest bardzo powolna, więc konwertujemy ją na minuty na cykl i wypisujemy tę wartość, która wynosi 17.14286.
Oto moje dane w formacie rozdzielanym tabulatorami, jeśli ktoś chce spróbować.
Dzięki za pomoc! Ten problem był trudny (dla mnie), ale świetnie się bawiłem!
Odpowiedzi:
Ciekawy projekt, który tam realizujesz! :-)
Z analizy sygnału POV jest to w rzeczywistości proste pytanie - i tak, masz rację, że użyłbyś FFT do tego problemu z estymacją częstotliwości.
Następnie, bardzo prosto, znajdź maksimum miejsca, w którym znajduje się twój PSD. Odcięta tego maksimum będzie odpowiadać twojej częstotliwości.
Zastrzegam Emptor, daję ci ogólne spojrzenie i podejrzewam, że wynikiem FFT w R będzie znormalizowana częstotliwość, w którym to przypadku musisz znać częstotliwość próbkowania (co robisz), aby ją przekonwertować na Hz. Jest wiele innych ważnych szczegółów, które pomijam, takich jak rozdzielczość częstotliwości, rozmiar FFT i fakt, że prawdopodobnie chcesz najpierw usunąć swój sygnał, ale dobrze będzie zobaczyć najpierw wykres.
EDYTOWAĆ:
Pozwól nam wziąć twój sygnał pod uwagę. Po tym, jak mam na myśli, wygląda to tak:
Możesz zobaczyć, jak to jest symetryczne. Jeśli zignorujesz ostatnią połowę i po prostu spojrzysz na pierwszą połowę i powiększysz, zobaczysz:
źródło
W przypadku przebiegu płynnego i stacjonarnego, zliczanie punktów próbnych między dodatnimi przejazdami o pewnej średniej wartości progowej da oszacowanie okresu. Spójrz na kilka okresów przekraczania progów, aby uzyskać bardziej średni szacunek lub wykryć dowolny trend.
źródło
Nie trzeba robić nic skomplikowanego: wystarczy zmierzyć czas trwania między pikami kształtu fali. To jest okres. Częstotliwość wynosi tylko 1 podzielone przez okres.
Przy około 8 cyklach w ciągu 2 godzin częstotliwość wynosi 4 cykle na godzinę lub około 1 mHz.
źródło