Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
# needs zoo for "rollapply"
# needs e1071 for Hamming window
# load data
th <- read.csv2('app_status_report_hourly.csv', header=F, sep=',', col.names=c('time','count'))
# normalize the data by the maximum
X <- th$count[1:length(th$count)]/max(th$count)
Xw <- rollapply(X,24,by=1,function(x){x*hamming.window(24)})
Fw <- apply(Xw,1,fft)
Fw_avg <- colMeans(t(abs(Fw)))
Fw_rem <- Fw - Fw_avg
Xw_rem <- apply(Xw,1,fft,inverse=T)
X_rem <- vector(mode=typeof(Xw_rem),length=length(X))
for(i in seq(1, dim(Xw_rem)[2])){
X_rem[i:(i+23)] <- X_rem[i:(i+23)] + Xw_rem[,i]
# plot(abs(X_rem),type='l') # the recovered with high frequency artifacts
plot(runmed(abs(X_rem), 25), type='l')
plot(runmed(abs(X), 25), type='l') #compare with simple median smoothing
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment