小波分析(by R)

小波分析

#采用WaveletComp包实现小波分析与画图，应注意小波分析数据需完整，存在缺失的话应先插值
library(WaveletComp)
data<-read.table('E:\\data\\data.txt',sep = "\t",encoding="UTF-8",header = T)
x = data$target
#数据从2020-01-01到2023-06-30，时间步长为一天
epoch.seq <- seq(from = as.POSIXct("2020-01-01"),
                 to = as.POSIXct("2023-06-30"), by = 3600*24)
#构建好小波分析的数据
my.data <- data.frame(date = epoch.seq, x = x)

#ticks与labels构建两个步长较大的序列用于小波分析画图
ticks <- seq(as.POSIXct("2020-01-01 00:00:00", format = "%F %T"),
             as.POSIXct("2023-06-30 23:00:00", format = "%F %T"), by = "month")
labels <- seq(as.Date("2020-01-01"), as.Date("2023-06-30"), by = "month")

#小波分析
my.w <- analyze.wavelet(my.data, "x",
                        loess.span = 0,#不对这个系列进行趋势分析
                        dt = 1, dj = 1/250,#每时间单位观察一次"x"(= 1/dt)。(这定义了时间单位。),dj决定沿周期轴的分辨率
                        lowerPeriod = 1,
                        upperPeriod = 1024,#定义小波变换中周期范围，以时间单位表示。只有在这个范围内的x周期才会被检测到。
                        #如果没有设置上周期和下周期，则默认值为lowerPeriod = 2*dt，upperPeriod = floor(nrow(my.data)/3)*dt
                        make.pval = TRUE, n.sim = 10)#通过10次模拟可以找到x中每个t的显著周期区域(图4中白线内的区域)。
                        
#设置默认路径，即后续图片的保存位置
setwd("E:\\output")
jpeg(filename =paste0("预测.jpg"),width=1000,height=1000,res=250)#res 是指 PPI（Pixel Per Inch） ，即图像的采样率（在图像中，每英寸所包含的像素数目），越大线条越粗
#画图
wt.image(my.w, periodlab = "periods (days)",label.time.axis = TRUE,
        color.key = "quantile", n.levels = 250,#n.levels是颜色级别的数量;默认值为100
        show.date = TRUE, date.format = "%F %T",
        spec.time.axis = list(at = ticks, labels = labels, las = 1),
        legend.params = list(lab = "wavelet power levels",mar = 4.7))#参数mar = 4.7控制颜色条到右边距的距离
dev.off()

jpeg(filename =paste0("显著性.jpg"),width=1000,height=1000,res=250)#res 是指 PPI（Pixel Per Inch） ，即图像的采样率（在图像中，每英寸所包含的像素数目），越大线条越粗

小波显著性

#显著性
wt.avg(my.w, siglvl = 0.01, sigcol = "red", sigpch = 20,maximum.level = 1.001*max(my.w$Power.avg),#各自的周期在1%的水平上是显著的。
       periodlab = "period (days)",legend.coords ="topright")
dev.off()

小波去噪

#去噪
jpeg(filename =paste0("溶解氧预测.jpg"),width=1000,height=1000,res=150)#res 是指 PPI（Pixel Per Inch） ，即图像的采样率（在图像中，每英寸所包含的像素数目），越大线条越粗

reconstruct(my.w,
            only.ridge = F, 
            lty = 1,
            lwd = c(1,3),
            col = c("grey","red"),
            show.date = TRUE,
            show.legend = T,
            legend.horiz=F,
            legend.text=c("原始曲线","去噪曲线"),
            date.format=NULL,
            date.tz=NULL,
            rescale = TRUE,
            main = "去噪图像",
            main.waves=NULL,
            main.rec=NULL,
            spec.time.axis = list(at = ticks, labels = labels, las = 1),
            legend.coords ="topright",
            ylim = c(min(data$target),max(data$target)))
dev.off()