Wavelet reconstruction of time series
I have been working on this topic currently, using the same paper. I show you code using an example dataset, detailing how I implemented the procedure of wavelet decomposition and reconstruction.
# Lets first write a function for Wavelet decomposition as in formula (1):
mo<-function(t,trans=0,omega=6,j=0){
dial<-2*2^(j*.125)
sqrt((1/dial))*pi^(-1/4)*exp(1i*omega*((t-trans)/dial))*exp(-((t-trans)/dial)^2/2)
}
# An example time series data:
y<-as.numeric(LakeHuron)
From my experience, for correct reconstruction you should do two things: first subject the mean to get a zero-mean dataset. I then increase the maximal scale. I mostly use 110 (although the formula in the Torrence and Compo suggests 71)
# subtract mean from data:
y.m<-mean(y)
y.madj<-y-y.m
# increase the scale:
J<-110
wt<-matrix(rep(NA,(length(y.madj))*(J+1)),ncol=(J+1))
# Wavelet decomposition:
for(j in 0:J){
for(k in 1:length(y.madj)){
wt[k,j+1]<-mo(t=1:(length(y.madj)),j=j,trans=k)%*%y.madj
}
}
#Extract the real part for the reconstruction:
wt.r<-Re(wt)
# Reconstruct as in formula (11):
dial<-2*2^(0:J*.125)
rec<-rep(NA,(length(y.madj)))
for(l in 1:(length(y.madj))){
rec[l]<-0.2144548*sum(wt.r[l,]/sqrt(dial))
}
rec<-rec+y.m
plot(y,type="l")
lines(rec,col=2)
As you can see in the plot, it looks like a perfect reconstruction:
wavelet package Inverse DWT fails to reconstruct series?
Following the help file ?dwt
, you can modify your script, such as:
library(wavelets)
set.seed(42)
x <- rnorm(100)
y <- idwt(dwt(x, n.levels=3, boundary="reflection", fast=FALSE))
plot(x, y)
abline(0,1)
Multilevel partial wavelet reconstruction with pyWavelets
I managed to write my own version of the wrcoef
function which appears to work as expected:
import pywt
import numpy as np
def wrcoef(X, coef_type, coeffs, wavename, level):
N = np.array(X).size
a, ds = coeffs[0], list(reversed(coeffs[1:]))
if coef_type =='a':
return pywt.upcoef('a', a, wavename, level=level)[:N]
elif coef_type == 'd':
return pywt.upcoef('d', ds[level-1], wavename, level=level)[:N]
else:
raise ValueError("Invalid coefficient type: {}".format(coef_type))
level = 4
X = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]
coeffs = pywt.wavedec(X, 'db1', level=level)
A4 = wrcoef(X, 'a', coeffs, 'db1', level)
D4 = wrcoef(X, 'd', coeffs, 'db1', level)
D3 = wrcoef(X, 'd', coeffs, 'db1', 3)
D2 = wrcoef(X, 'd', coeffs, 'db1', 2)
D1 = wrcoef(X, 'd', coeffs, 'db1', 1)
print A4 + D4 + D3 + D2 + D1
# Results:
[ 9.99200722e-16 1.00000000e+00 2.00000000e+00 3.00000000e+00
4.00000000e+00 5.00000000e+00 6.00000000e+00 7.00000000e+00
8.00000000e+00 9.00000000e+00 1.00000000e+01 1.10000000e+01
1.20000000e+01 1.30000000e+01 1.40000000e+01 1.50000000e+01
1.60000000e+01 1.70000000e+01]
Related Topics
Navlistpanel: Make Tabs Sequentially Active in Shiny App
Solving a System of Nonlinear Equations in R
R - Reading Lines from a .Txt-File After a Specific Line
R Cumulative Sum with a Condition and a Reset
How to Rename Element's List Indexed by a Loop in R
How to Add Random 'Na's into a Data Frame
R Dplyr Join on Range of Dates
Ggplot Boxplot - Length of Whiskers with Logarithmic Axis
Generating a Color Legend with Shifted Labels Using Ggplot2
Change Color Median Line Ggplot Geom_Boxplot()
Extracting Zip+CSV File from Attachment W/ Image in Body of Email
R Histogram with Multiple Populations
Select Last Row by Group for All Columns Data.Table
Gap in Polar Time Plot - How to Connect Start and End Points in Geom_Line or Remove Blank Space
Adjusting the Width of Legend for Continuous Variable
R Obtaining Rownames Date Using Quantmod
Tidyverse Not Loaded, It Says "Namespace 'Vctrs' 0.2.0 Is Already Loaded, But >= 0.2.1 Is Required"