Określanie częstotliwości i okresu fali

11

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

Moja fala

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

znormalizowany wykres

Na razie w porządku. Oto wykres gęstości widmowej:

gęstość widmowa

Następnie powiększamy po lewej stronie wykresu i zaznaczamy najwyższy wskaźnik (czyli 70) zieloną linią:

powiększyć wykres spektralny

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!

Aaron Patterson
źródło
Aaron, myślę, że najlepszą rzeczą jest, aby umieścić link do pliku danych (jako tekst lub coś takiego) na Dropbox, aby móc go pobrać i dać odpowiedź. W przeciwnym razie będzie dużo do przodu i do tyłu. Nie widzę liczb na skrajnym lewym końcu. :-) (Podaj także częstotliwość próbkowania - czyli jak często odczytujesz temperaturę).
Spacey
O przepraszam. Dane zawierają temperaturę w stopniach C, przeliczam na stopnie F dla wykresu. Są to jednak poprawne dane (kolumna „temp”).
Aaron Patterson
Problem z pomiarem częstotliwości w ten sposób polega na tym, że jeśli uzyskasz znaczną zmienność między cyklami, wówczas znacznie trudniej będzie ustalić średnią częstotliwość - szczyty będą się rozmazywać - podczas gdy samo zliczanie czasu między wycieczkami pozwoli ci dobrze uśrednić rzeczy (a także oblicza std dev itp.). Zastosowanie metody FFT byłoby bardziej uzasadnione, gdyby było dużo hałasu, ale wydaje się, że tak nie jest.
Daniel R Hicks
+1 za publikowanie, kod, dane, wykresy i link do github.
nibot
@DanielRHicks W tym konkretnym przypadku nie sądzę, żeby to miało znaczenie, ale tak, FFT da ci średnią z nich wszystkich, podczas gdy gdybyśmy zrobili coś w rodzaju przejścia przez zero, mierzylibyśmy czas trwania każdego cyklu (częstotliwość), i wtedy możemy ustalić, czy chcemy wziąć średnią, medianę, tryb itp. Dobra uwaga!
Spacey,

Odpowiedzi:

7

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.

real2+imag2

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:

wprowadź opis zdjęcia tutaj

x[n]

Nfft=103600=36000.

fs=0.5Hz

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

wprowadź opis zdjęcia tutaj

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:

wprowadź opis zdjęcia tutaj

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

1(9.7222e004)60=17.14

Spacey
źródło
@AaronPatterson Zredagowałem post, proszę zobaczyć. Możesz także dodawać zdjęcia bezpośrednio do oryginalnego postu. :-). Dodaj obraz otrzymanego wyniku PSD.
Spacey,
1
Nie do końca poprawne, jeśli częstotliwość okaże się między przedziałami wyników FFT.
hotpaw2
@ hotpaw2 Dlatego ostrzegłem OP, że daję ogólne perspektywy i dlatego muszę zobaczyć fabułę. Mimo to edytowałem, aby dodać dodatkowe zastrzeżenia.
Spacey
1
@AaronPatterson Bez problemu, chętnie pomogę. Jeśli chodzi o książki, spójrz na „Zrozumienie DSP” Richarda Lyonsa - jest to szybka książka na początek.
Spacey
1
1.3x105
4

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.

hotpaw2
źródło
3

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.

nibot
źródło
3
Jak mogę to zrobić programowo?
Aaron Patterson