R语言利用wmtsa进行MODWT分析的一些问题

这篇文章的作者利用R语言的wmtsa包进行了MODWT分析 ¹，我测试了一下发现有点问题。关于MODWT网上已经有很多帖子介绍了，这里就不说明。

安装并载入wmtsa包

安装命令：install.packages(‘wmtsa’)
载入命令：library(wmtsa)

wavMODWT

这是执行maximal overlap discrete wavelet transform的函数，help(wavMODWT)可以看到如下信息：

用法

wavMODWT(x, wavelet="s8", n.levels=ilogb(length(x), base=2),
         position=list(from=1,by=1,units=character()), units=character(),
         title.data=character(), documentation=character(), keep.series=FALSE)

重要参数

x：输入的实数时间序列
keep.series：逻辑值，‘TRUE’或者‘FALSE’，决定是否在输出对象中保留原始时间序列，默认FALSE
n.levels：分解层数，默认为最大可分解层数即$\lfloor {\log }_{2}n\rfloor$，其中$n$为数据长度
wavelet：分解采用的小波函数，默认为s8

运行示例

linchirp <- make.signal("linchirp", n=1024)
result   <- wavMODWT(linchirp, wavelet="s8", n.levels=5, keep.series=TRUE)
plot(wavShift(result))
eda.plot(result)

代码的大致意思是生成长度1024的线性调频信号，分解5层，并且保留原始时间序列

这是summary的绘制，包括系数、能量分布和分位数

reconstruct

该函数只能一次性重构所有信号而非多尺度分解，这里用不上。另外，MODWT因为放弃了正交特性，直接逆变换回去不是多尺度分析想要的结果，因而需要特别的函数。

wavMRD

该函数可以实现多尺度分解，如果没有特殊要求，用法极为简单：rec=wavMRD(result)。rec是一个对象，包括所有分解层数；如果要获取某一层对象，引用方式为一对中括号rec[‘D1’]或者rec[1]，其中引号内为目标层；如果要获取某一层的数据，引用方式为两对中括号rec[[‘D1’]]或者rec[[1]]。通过matrix_rec=as.matrix(rec)可以将对象转换为二维矩阵更便于使用。对于该对象可以进行的操作还有：boxplot绘制箱图；crystal.names获取尺度名称；plot直接绘制数据。
遗憾的是wavMRD分解效果在实际数据中并不好，下图是本人的一个例子，逐渐分解下去，可以看到边界明显的端点效应，以及大尺度下的分解结果相比原始数据拟合效果非常差，这是因为wavMRD里面没有reflect参数，无法进行reflect延拓。

hwt=wavMODWT(h)
rech=wavMRD(hwt)
mah=as.matrix(rech)
m=dim(mah)[1];n=dim(mah)[2];
reci=h
par(mfrow=c(3,4))
for(i in (1:n))
{
	plot(t,h,type='o')	
	points(t,reci,col='red',type='l',lwd=4);
	#lines(t,reci,col='red')
	reci=reci-mah[,i]
}

wavMRDSum

wavMRDSum虽然是累加多个尺度分解结果的函数，但是通过设置目标尺度为一个即可获取单层分解结果。

用法

wavMRDSum(x, wavelet="s8", levels=1, xform="modwt", reflect=TRUE,
         keep.smooth=TRUE, keep.details=TRUE)

重要参数

x：输入的实数时间序列
keep.details：是否保留细节信号，默认为TRUE
keep.smooth：是否保留光滑信号，默认为TRUE
levels：目标尺度，默认为1；在不添加反射延拓条件时的最大层数为：$\lfloor {\log }_{2}n\rfloor$，其中$n$为数据长度；在添加反射延拓条件时的最大层数为：$\lfloor {\log }_2\frac{n-1}{L-1}\rfloor$，其中$L$为小波滤波器长度。
reflect：是否进行反射延拓，默认为TRUE
wavelet：小波函数
xform：分析方法，可选DWT和MODWT，默认MODWT

运行示例

par(mfrow=c(3,4))
for(i in (1:11))
{
	plot(t,h,type='o')
	points(t,wavMRDSum(h, levels=i,wavelet='s2'),col='red',type='l',lwd=4);
}

结论

wavMRD可以进行多尺度分解，但是没有反射延拓参数，所以分解结果较差。wavMRDSum有反射延拓参数但是不知道为什么对于分解层数的限制和wavMRD不一样，也就是说和系数能具有的最大层数不一致？个人觉得这个包的设计哲学是有问题的，如果要通过wavMRDSum来实现反射延拓的多尺度分析，那么就要进行多次小波变换多次重构，而不是一次小波变换多次重构，运算效率降低。我认为解决方案是wavMRD中添加反射延拓参数，或者将wavMRDSum的输入变为小波变换的结果，这样可以减少运算量。最后，我还是推荐利用MATLAB进行MODWT多尺度分析，它对于反射延拓下的分解层数没有限制，和系数层数一致，有利于编程进行分析。

---

1. Watson P J . Identifying the Best Performing Time Series Analytics for Sea Level Research[M]// Time Series Analysis and Forecasting. Springer International Publishing, 2016. ↩︎