1

我需要一些帮助来进行时间序列分析，特别是带有汉明窗平滑的快速傅里叶变换。

TL;博士

  1. fftper() 是 R 中带有汉明窗平滑的 FFT 的适当函数吗？

  2. 如何提取或生成 fftper 输出图的频率值？

  3. 我正在寻找每日周期，因此要将频率数据转换回“时间”变量，我是否将频率值除以 1/24？（例如这里的麻疹示例：http ://web.stanford.edu/class/earthsys214/notes/series.html ）

长版

我已经为一堆个体动物的声学检测数据加上了时间戳。对于每个人，我将检测次数按小时分组，并使用频率为 24 的 ts() 将其转换为时间序列（查看检测中的每日模式）。有了这些数据，我想应用汉明窗平滑，然后进行 FFT 并生成这些数据的周期图。我还想提取频率值（x 轴）并将它们从频率转换为时间段（小时）。

我已经设法在 FFT 上生成周期图，并fftper()使用TSSS package. 自动生成的图看起来不错。我现在想提取图中使用的频率值（x 轴值）和功率值（y 轴值），因此我可以将频率（x 轴值）数据转换回时间变量（即使用我认为频率/(1/24)?)，然后用ggplot. fftper生成一个spg结构如下的对象：

List of 4
 $ period         : num [1:65] 10.43 1.95 2.36 2.57 1.9 ...
 $ smoothed.period: num [1:65] 0.815 0.601 0.364 0.374 0.348 ...
 $ log.scale      : chr "TRUE"
 $ tsname         : chr "hourly_ts"
 - attr(*, "class")= chr "spg"

我可以提取 y 轴值（平滑周期值或功率），FFTpower <- FFT[["smoothed.period"]]但我看不到 x 轴值的存储位置或弄清楚如何生成它们。

有任何想法吗？提前致谢！

虚拟数据：

#Data
df <- read.table(text =
               "timestampUTC    ID
'2017-10-02 19:23:27'    47280
'2017-10-02 19:26:48'    47280
'2017-10-02 19:27:23'    47280
'2017-10-02 19:31:46'    47280
'2017-10-02 23:52:15'    47280
'2017-10-02 23:53:26'    47280
'2017-10-02 23:55:13'    47280
'2017-10-03 19:53:50'    47280
'2017-10-03 19:55:23'    47280
'2017-10-03 19:58:26'    47280
'2017-10-04 13:15:13'    47280
'2017-10-04 13:16:42'    47280
'2017-10-04 13:21:39'    47280
'2017-10-04 19:34:54'    47280
'2017-10-04 19:55:28'    47280
'2017-10-04 20:08:23'    47280
'2017-10-04 20:21:43'    47280
'2017-10-05 04:55:48'    47280
'2017-10-05 04:57:04'    47280
'2017-10-05 05:18:40'    47280
'2017-10-07 21:24:19'    47280
'2017-10-07 21:25:36'    47280
'2017-10-07 21:29:25'    47280", header = T)

代码：

#convert datetime
df$timestampUTC<-as.POSIXct(df$timestampUTC, format = "%Y-%m-%d %H:%M:%S", tz="UTC")

#keep only datetime column and add second column with frequency of 1
df<-df %>%
  select(timestampUTC)
df<-data.frame(df,Frequency=1)

#bin into hours
hourly_detections <- df %>%
  mutate(processed_hour = floor_date(timestampUTC, "hour")) %>%
  group_by(processed_hour)%>%
  summarise(count = sum(Frequency))

#set time frame using max and min hours
time_frame <- as_datetime(c(min(floor_date(df$timestampUTC,"hour")),(max(ceiling_date(df$timestampUTC,"hour"))-1)),tz="Australia/Sydney")

#combine detection hour and non-detections hour dfs
all_hours <- data.frame(processed_hour = seq(time_frame[1], time_frame[2], by = "hour"))

#build df with every hour and set count to 0 for 'new' hours
hourly_detections <- hourly_detections %>%
  right_join(all_hours, by = "processed_hour") %>%
  mutate(count = ifelse(test = is.na(count),yes  = 0,no   = count))
hourly_detections<-hourly_detections[order(hourly_detections$processed_hour),]

#set up time series
hourly_ts <- ts(hourly_detections$count, start= min(hourly_detections$processed_hour), frequency=24) 

#FFT with hamming widow smoothing
FFT<-fftper(hourly_ts, window = 2, plot = TRUE)

#extract y (power) values
FFTPower<-FFT[["smoothed.period"]]

#extract x values?
4

0 回答 0