Określanie częstotliwości i okresu fali

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!

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.

$real^2 + imag^2$

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:

$x[n]$

$N_{fft} = 10 * 3600 = 36000.$

$f_s = 0.5 Hz$

$x[n]$$X(f)$$10log_{10}(|X(f)|^{2})$

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:

$\frac{F_s}{N_{fft}} = 1.3889e-005 Hz$$0 * 1.3889e-005 = 0 Hz$$1 * 1.3889e-005 = 1.3889e-005 Hz$$70 * 1.3889e-005 = 9.7222e-004 Hz$

$\frac{1}{(9.7222e-004)*60} = 17.14$

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.

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.