I need some help doing time series analysis, specifically Fast Fourier transformations with hamming window smoothing.
TL;DR
is fftper() the appropriate function for FFT with hamming window smoothing in R?
How can I extract or generate the frequency values of the fftper output plot?
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?