Я собираю данные о температуре из холодильника. Данные выглядят как волна. Я хотел бы определить период и частоту волны (чтобы я мог измерить, если модификации холодильника имеют какой-либо эффект).
Я использую R, и я думаю, что мне нужно использовать FFT для данных, но я не уверен, куда идти дальше. Я очень плохо знаком с R и анализом сигналов, поэтому любая помощь будет принята с благодарностью!
Вот волна, которую я создаю:
Вот мой код R до сих пор:
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))
Я разместил код R вместе с базой данных SQLite здесь .
Вот график нормализованного графа (с удаленным средним):
Все идет нормально. Вот график спектральной плотности:
Затем мы увеличиваем левую часть графика и отмечаем самый высокий индекс (который равен 70) зеленой линией:
Наконец мы рассчитываем частоту волны. Эта волна очень медленная, поэтому мы конвертируем ее в минуты за цикл и выводим значение, которое равно 17.14286.
Вот мои данные в формате с разделителями табуляции, если кто-то еще хочет попробовать.
Спасибо за помощь! Эта проблема была трудной (для меня), но я прекрасно провел время!
Ответы:
Интересный проект у вас там происходит! :-)
С точки зрения анализа сигналов это на самом деле простой вопрос - и да, вы правы, что вы использовали бы БПФ для этой задачи оценки частоты.
Тогда, очень просто, найдите максимум того, где находится ваш PSD. Абсцисса этого максимума будет соответствовать вашей частоте.
Будьте бдительны, я даю вам общий обзор, и я подозреваю, что результатом БПФ в R будет нормализованная частота, и в этом случае вам придется знать частоту дискретизации (что вы делаете), чтобы преобразовать ее обратно. в Гц. Есть много других важных деталей, которые я пропускаю, таких как разрешение по частоте, размер БПФ, и тот факт, что вы, вероятно, хотите сначала деинформировать свой сигнал, но сначала будет полезно увидеть график.
РЕДАКТИРОВАТЬ:
Позвольте нам принять ваш сигнал во внимание. После того, как я обманываю это, это выглядит так:
Вы можете видеть, как это симметрично. Если вы игнорируете последнюю половину и просто смотрите на первую половину и увеличиваете масштаб, вы можете увидеть это:
источник
Для сигнала, такого плавного и стационарного, подсчет точек выборки между положительными проходящими пересечениями некоторого среднего порогового значения даст оценку периода. Посмотрите на несколько периодов пересечения порогов, чтобы получить более среднюю оценку или обнаружить любую тенденцию.
источник
Нет необходимости делать что-либо сложное: просто измерьте длительность между пиками формы сигнала. Это период. Частота всего 1, деленная на период.
При примерно 8 циклах в течение 2 часов частота составляет 4 цикла в час или около 1 МГц.
источник