对R中的时间序列执行傅立叶分析 [英] Perform Fourier Analysis to a Time Series in R
问题描述
我想使用R对时间序列执行傅立叶变换.我想:
I would like to perform fourier transform to a time series using R. I would like to:
- 获取5至18次谐波的总和
- 绘制每个波浪
- 并输出为csv文件.
以下是数据的链接: 链接到数据
Here's the link to the data: Link to Data
这是我的初始代码.
dat <- read.csv("Baguio.csv",header=FALSE)
y <- dat$V1
ssp <-spectrum(y)
t <- 1:73
per <- 1/ssp$freq[ssp$spec==max(ssp$spec)]
reslm <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t))
rg <- diff(range(y))
#blue dashed line
plot(y~t,ylim=c(min(y)-0.1*rg,max(y)+0.1*rg))
lines(fitted(reslm)~t,col=4,lty=2)
#green line 2nd harmonics
reslm2 <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t))
lines(fitted(reslm2)~t,col=3)
有没有一种方法可以简化此代码?如果必须达到18次谐波，该方程将变得非常长.另外，我仍然不知道如何在此处添加谐波.
Is there a way to simplify this code? If I have to get to the 18th harmonics the equation becomes very very long. Also, I still do not know how to add the harmonics here.
在此先感谢
推荐答案
更简单的解决方案是使用快速傅立叶变换(fft)
A much easier solution is to use the Fast Fourier Transform (fft)
dat <- read.csv("Baguio.csv", header=FALSE)
y <- dat$V1
t <- 1:73
rg <- diff(range(y))
nff = function(x = NULL, n = NULL, up = 10L, plot = TRUE, add = FALSE, main = NULL, ...){
#The direct transformation
#The first frequency is DC, the rest are duplicated
dff = fft(x)
#The time
t = seq(from = 1, to = length(x))
#Upsampled time
nt = seq(from = 1, to = length(x)+1-1/up, by = 1/up)
#New spectrum
ndff = array(data = 0, dim = c(length(nt), 1L))
ndff[1] = dff[1] #Always, it's the DC component
if(n != 0){
ndff[2:(n+1)] = dff[2:(n+1)] #The positive frequencies always come first
#The negative ones are trickier
ndff[length(ndff):(length(ndff) - n + 1)] = dff[length(x):(length(x) - n + 1)]
}
#The inverses
indff = fft(ndff/73, inverse = TRUE)
idff = fft(dff/73, inverse = TRUE)
if(plot){
if(!add){
plot(x = t, y = x, pch = 16L, xlab = "Time", ylab = "Measurement",
main = ifelse(is.null(main), paste(n, "harmonics"), main))
lines(y = Mod(idff), x = t, col = adjustcolor(1L, alpha = 0.5))
}
lines(y = Mod(indff), x = nt, ...)
}
ret = data.frame(time = nt, y = Mod(indff))
return(ret)
}
然后，我们需要调用res，将时间序列传递为x，将谐波数传递为n，并将上采样(因此我们将时间点绘制在原始采样点旁边)作为up.
Then we need to call res, passing the timeseries as x, the number of harmonics as n and the upsampling (so we plot points in time beside the original ones) as up.
png("res_18.png")
res = nff(x = y, n = 18L, up = 100L, col = 2L)
dev.off()
要获得5至18次谐波之和，只是序列之间的差异
To get the sum of the 5th to the 18th harmonics it's simply a difference between series
sum5to18 = nff(x = y, n = 18L, up = 10L, plot = FALSE)
sum5to18$y = sum5to18$y - nff(x = y, n = 4L, up = 10L, plot = FALSE)$y
png("sum5to18.png")
plot(sum5to18, pch = 16L, xlab = "Time", ylab = "Measurement", main = "5th to 18th harmonics sum", type = "l", col = 2)
dev.off()
添加参数add和col允许我们绘制具有特定颜色的多条波浪
Adding the arguments add and col allow us to plot multiple waves as well, with specific colors
colors = rainbow(36L, alpha = 0.3)
nff(x = y, n = 36L, up = 100L, col = colors[1])
png("all_waves.png")
for(i in 1:18){
ad = ifelse(i == 1, FALSE, TRUE)
nff(x = y, n = i, up = 100L, col = colors[i], add = ad, main = "All waves up to 18th harmonic")
}
dev.off()
有没有办法提取每个系列的数据，然后另存为csv文件.因此，在此示例中，我应该为18个波形设置18个csv文件.
Is there a way so extract the data of each series then save as a csv file. So in this example, I should have 18 csv files for the 18 waves.
我编辑了代码，以允许谐波为0(基本上是平均值)，因此现在您将单独的波提取为:
I edited the code to allow a 0 harmonic (basically a mean), so now you extract the separate waves as:
sep = array(data = NA_real_, dim = c(7300L, 2 + 18), dimnames = list(NULL, c("t", paste0("H", 0:18))))
sep[,1:2] = as.matrix(nff(x = y, n = 0, up = 100L, plot = FALSE))
for(i in 1:18L){
sep[,i+2] = nff(x = y, n = i, up = 100L, plot = FALSE)$y - nff(x = y, n = i-1, up = 100L, plot = FALSE)$y
}
然后，您可以使用write.table编写一个csv文件.
Then you can use write.table to write a csv file.
这篇关于对R中的时间序列执行傅立叶分析的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
