3

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

为了重构信号，我使用了Torrence & Compo 经典文本中的公式（11） ，但结果与原始信号无关。我特别关心小波变换的实部与尺度之间的划分，这一步完全扭曲了结果。另一方面，如果我只是对所有尺度的实部求和，结果与原始时间序列非常相似，但值稍宽（原始序列范围~ [-0.2, 0.5]，重构序列范围〜[-0.4,0.7]）。

我想知道是否有人能说出一些实用的程序、公式或算法来重建原始时间序列。我已经阅读了 Torrence 和 Compo (1998)、Farge (1992) 和其他书籍的论文，它们都有不同的公式，但没有人真正帮助我。

4

1 回答 1

6

我目前一直在研究这个主题，使用同一篇论文。我向您展示了使用示例数据集的代码，详细说明了我如何实现小波分解和重建的过程。

# 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)

根据我的经验，为了正确重建，您应该做两件事：首先对均值进行处理以获得零均值数据集。然后我增加最大比例。我主要使用 110（尽管 Torrence 和 Compo 中的公式建议使用 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)

正如您在情节中看到的那样，它看起来像是一个完美的重建：

原始级数和小波重建

于 2014-03-02T18:17:31.700 回答