确定波的频率和周期

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

#1 楼

您正在进行的有趣项目！ :-)

从信号分析POV来看，这实际上是一个简单的问题-是的，您是对的，您可以将FFT用于此频率估计问题。

我对R并不熟悉，但是您实际上想要做的是对温度信号进行FFT。由于您的信号是真实的，因此您将获得复杂的结果作为FFT的输出。现在，您可以采用复杂FFT结果的绝对幅度平方（因为我们不在乎相位）。即$ real ^ 2 + imag ^ 2 $。这是您的功率谱密度。

然后，非常简单地找到PSD所在的最大值。该最大值的横坐标将对应于您的频率。

Caveat Emptor，我给您一个大致的看法，我怀疑R中FFT的结果将是归一化频率，在这种情况下，您将必须知道您的采样率，（ ），以便将其转换回Hz。我遗漏了许多其他重要的细节，例如您的频率分辨率，FFT大小，以及您可能想先去除信号平均值的事实，但是最好先查看图表。

编辑：

让我们考虑您的信号。在我消除均值之后，它看起来像这样：



您要消除均值，因为技术上DC偏置的频率为0 Hz，您可以不想让它淹没其他感兴趣的频率。让我们称此为非均值信号$ x [n] $。

现在，当您进行FFT运算时，必须注意一些细节。我使用的FFT长度是信号长度的10倍，因此我们可以说$ N_ {fft} = 10 * 3600 = 36000. $任何设置为大于原始信号大小的FFT大小都会对您的信号进行内插FFT产生频域。尽管这不会从信息论POV中为您添加任何信息，但它确实可以帮助您更好地识别出真正的峰值位于何处，尤其是在两个频点之间的情况下。为了获得真正更好的频率分辨率，您想获取更多数据。 （即，让传感器运行更长的时间）。

从文本文件来看，我发现您的传感器每2秒获取一次温度读数。这意味着您的采样率，即$ f_s = 0.5 Hz该结果将是复杂的，并且实际上是共轭对称的，但这在这里并不重要。 （这只是出于您的目的，其中一半是多余的）。现在，如果我们取$ 10log_ {10}（| X（f）| ^ {2}）$，这就是功率谱密度（以对数标度表示）。结果如下所示：



您可以看到它是如何对称的。如果您忽略后半部分，而只是看一下前半部分并进行放大，则会看到以下内容：但是实际频率指数是多少？这是您想要的频率分辨率的地方。为了进行计算，我们只需取\\ frac {F_s} {N_ {fft}} = 1.3889e-005 Hz $。这仅表示每个索引代表的频率是该数字的倍数。索引0为$ 0 * 1.3889e-005 = 0 Hz $。索引1是$ 1 * 1.3889e-005 = 1.3889e-005 Hz $。索引70为$ 70 * 1.3889e-005 = 9.7222e-004 Hz $。

我们快完成了。现在您已经有了这个数字，我们可以将其转换为更可口的东西。这个数字只是说您的信号每秒完成9.7222e-004个周期。因此，我们可以问一下，计算一个周期需要多少分钟？因此，每个周期$ \ frac {1} {（9.7222e-004）* 60} = 17.14 $分钟。

评论

---

$ \ begingroup $
@AaronPatterson我已经编辑了帖子，请参阅。另外，您可以将图片直接添加到原始帖子中。 :-)。请添加您获得的PSD结果的图像。
$ \ endgroup $
–太空
2012年6月2日19:18

---

$ \ begingroup $
如果频率结果在FFT结果块之间，则不完全正确。
$ \ endgroup $
– hotpaw2
2012年6月2日21:15

---

$ \ begingroup $
@ hotpaw2这就是为什么我警告OP我要给出总体观点，以及为什么我需要看情节的原因。都一样，我进行了编辑以添加其他警告。
$ \ endgroup $
–太空
2012年6月2日21:35

---

$ \ begingroup $
@AaronPatterson没问题，很高兴为您提供帮助。至于书籍，请看Richard Lyons的“ Understanding DSP”-这是一本入门的快速书。
$ \ endgroup $
–太空
2012年6月4日20:57

---

$ \ begingroup $
@Frank这只是一个指数。例如，$ 1.3 \ text {x} 10 ^ {-5} $。
$ \ endgroup $
–太空
2014年1月20日12:54

---