3

I'm trying to reconstruct the original time series from a Morlet's wavelet transform. I'm working in R, package Rwave, function cwt. The result of this function is a matrix of n*m (n=period, m=time) containing complex values.

To reconstruct the signal I used the formula (11) in Torrence & Compo classic text, but the result has nothing to do with the original signal. I'm specially concerned with the division between the real part of the wavelet transform and the scale, this step distorts completely the result. On the other hand, if I just sum the real parts over all the scales, the result is quite similar to the original time series, but with slightly wider values (the original series ranges~ [-0.2, 0.5], the reconstructed series ranges ~ [-0.4,0.7]).

I'm wondering if someone could tell of some practical procedure, formula or algorithm to reconstruct the original time series. I've already read the papers of Torrence and Compo (1998), Farge (1992) and other books, all with different formulas, but no one really help me.

Konrad Rudolph
  • 530,221
  • 131
  • 937
  • 1,214
user3369539
  • 35
  • 1
  • 4

1 Answers1

7

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:

original series and wavelet reconstruction

ClimateUnboxed
  • 7,106
  • 3
  • 41
  • 86
DatamineR
  • 10,428
  • 3
  • 25
  • 45
  • Many thanks Rstudent, the code works fine. I am still wondering why the cwt function in Rwave doesn't produce the same results, but your code is enough for me now. – user3369539 Mar 04 '14 at 18:02
  • 1
    this answer is one of the most useful I have found on stackoverflow. I find most packages on R/matlab etc do not provide wavelet/inverse wavelet transforms that combine to the identity function. Thank you! – ClimateUnboxed Jan 31 '18 at 09:23
  • Hello @DatamineR I see in Torrence that the constant 0.2144548 is not trivial to compute. Could you share some code to use with other wt's, or maybe you know where to find them? Thank you in advance – fina Mar 04 '23 at 21:12