FFT with hamming window smoothing and extracting data from automatically generated plot

417 views Asked by At

I need some help doing time series analysis, specifically Fast Fourier transformations with hamming window smoothing.

TL;DR

  1. is fftper() the appropriate function for FFT with hamming window smoothing in R?

  2. How can I extract or generate the frequency values of the fftper output plot?

  3. I am looking for daily cycles, so to transform the frequency data back into a 'time' variable, do I divide the frequency values by 1/24? (e.g. measles example here: http://web.stanford.edu/class/earthsys214/notes/series.html)

Long version

I have timestamped acoustic detection data for a bunch of individual animals. For each individual, I have binned the number of detections by hour and converted this to a time series using ts() with a frequency of 24 (looking at daily patterns in detections). With these data, I want to apply hamming window smoothing and then a FFT and generate a periodogram of these data. I also want extract the frequency values (x axis) and convert these from frequency to time period (hours).

I have managed to generate the periodogram on the FFT and hamming windowed data this using the fftper() function in the TSSS package. The automatically generated plots look right. I now want to extract the frequency values (x axis values) and the power values (y axis values) used in the plot, so I can transform the frequency (x axis values) data back into a time variable (i.e. using I think frequency/(1/24)?), and then plot it nicely with ggplot. fftper generates an spg object that is structured like this:

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"

I can extract the y-axis values (smoothed period values or power) with FFTpower <- FFT[["smoothed.period"]] but I can't see where the x-axis values are stored or figure out how to generate them.

Any ideas? Thanks in advance!

Dummy data:

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

Code:

#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?
0

There are 0 answers