时间序列的小波重构 [英] Wavelet reconstruction of time series

问题描述

我正在尝试从 Morlet 的小波变换重建原始时间序列.我在 R 中工作，包 Rwave，函数 cwt.该函数的结果是一个包含复数值的 n*m (n=period, m=time) 矩阵.

为了重建信号，我使用了

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.

解决方案

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:

这篇关于时间序列的小波重构的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！