如何扩展pyWavelets以处理N维数据?
这可能是个适合其他论坛的问题,如果是这样,请告诉我。我注意到只有14个人关注小波这个标签。
我这里有一种优雅的方法,可以在pywt(pyWavelets包)中将小波分解扩展到多个维度。如果你安装了pywt,这个代码应该可以直接运行。测试1展示了一个3D数组的分解和重组。只需要增加维度的数量,代码就能处理4、6甚至18维的数据。
我在这里替换了pywt.wavedec和pywt.waverec这两个函数。在fn_dec中,我展示了新的wavedec函数是如何像旧的那样工作的。
不过有一个问题:它将小波系数表示为与数据形状相同的数组。因此,根据我对小波的有限了解,我只能用它处理Haar小波。像DB4这样的其他小波在这个严格的边界上会出现系数溢出(但目前的系数表示方法是一个数组列表[CA, CD1 ... CDN],这不是问题)。另一个问题是,我只在2^N边缘的立方体数据上使用过这个方法。
从理论上讲,我认为应该可以确保不发生“溢出”。关于这种小波分解和重组的算法在《C语言数值 Recipes》中有讨论,作者是William Press、Saul A Teukolsky、William T. Vetterling和Brian P. Flannery(第二版)。虽然这个算法假设在边缘进行反射,而不是其他形式的边缘扩展(比如零填充),但这个方法足够通用,可以适用于其他扩展形式。
有没有人建议如何将这个工作扩展到其他小波上?
注意:这个问题也发布在 http://groups.google.com/group/pywavelets
谢谢,
Ajo
import pywt
import sys
import numpy as np
def waveFn(wavelet):
if not isinstance(wavelet, pywt.Wavelet):
return pywt.Wavelet(wavelet)
else:
return wavelet
# given a single dimensional array ... returns the coefficients.
def wavedec(data, wavelet, mode='sym'):
wavelet = waveFn(wavelet)
dLen = len(data)
coeffs = np.zeros_like(data)
level = pywt.dwt_max_level(dLen, wavelet.dec_len)
a = data
end_idx = dLen
for idx in xrange(level):
a, d = pywt.dwt(a, wavelet, mode)
begin_idx = end_idx/2
coeffs[begin_idx:end_idx] = d
end_idx = begin_idx
coeffs[:end_idx] = a
return coeffs
def waverec(data, wavelet, mode='sym'):
wavelet = waveFn(wavelet)
dLen = len(data)
level = pywt.dwt_max_level(dLen, wavelet.dec_len)
end_idx = 1
a = data[:end_idx] # approximation ... also the original data
d = data[end_idx:end_idx*2]
for idx in xrange(level):
a = pywt.idwt(a, d, wavelet, mode)
end_idx *= 2
d = data[end_idx:end_idx*2]
return a
def fn_dec(arr):
return np.array(map(lambda row: reduce(lambda x,y : np.hstack((x,y)), pywt.wavedec(row, 'haar', 'zpd')), arr))
# return np.array(map(lambda row: row*2, arr))
if __name__ == '__main__':
test = 1
np.random.seed(10)
wavelet = waveFn('haar')
if test==0:
# SIngle dimensional test.
a = np.random.randn(1,8)
print "original values A"
print a
print "decomposition of A by method in pywt"
print fn_dec(a)
print " decomposition of A by my method"
coeffs = wavedec(a[0], 'haar', 'zpd')
print coeffs
print "recomposition of A by my method"
print waverec(coeffs, 'haar', 'zpd')
sys.exit()
if test==1:
a = np.random.randn(4,4,4)
# 2 D test
print "original value of A"
print a
# decompose the signal into wavelet coefficients.
dimensions = a.shape
for dim in dimensions:
a = np.rollaxis(a, 0, a.ndim)
ndim = a.shape
#a = fn_dec(a.reshape(-1, dim))
a = np.array(map(lambda row: wavedec(row, wavelet), a.reshape(-1, dim)))
a = a.reshape(ndim)
print " decomposition of signal into coefficients"
print a
# re-composition of the coefficients into original signal
for dim in dimensions:
a = np.rollaxis(a, 0, a.ndim)
ndim = a.shape
a = np.array(map(lambda row: waverec(row, wavelet), a.reshape(-1, dim)))
a = a.reshape(ndim)
print "recomposition of coefficients to signal"
print a
1 个回答
首先,我想给你推荐一个已经实现了单层多维变换的函数(来源)。这个函数会返回一个包含n维系数数组的字典。系数是通过描述变换类型(近似/细节)的键来访问的,这些变换应用于每个维度。
举个例子,对于二维的情况,结果是一个包含近似和细节系数数组的字典:
>>> pywt.dwtn([[1,2,3,4],[3,4,5,6],[5,6,7,8],[7,8,9,10]], 'db1')
{'aa': [[5.0, 9.0], [13.0, 17.0]],
'ad': [[-1.0, -1.0], [-1.0, -1.0]],
'da': [[-2.0, -2.0], [-2.0, -2.0]],
'dd': [[0.0, 0.0], [0.0, -0.0]]}
其中,aa
是对两个维度都应用了近似变换的系数数组(LL),而da
是对第一个维度应用了细节变换,对第二个维度应用了近似变换的系数数组(HL)(可以和dwt2的输出进行比较)。
基于这个,扩展到多层的情况应该是相对简单的。
这是我对分解部分的看法:https://gist.github.com/934166.
我还想提到你在问题中提到的一个问题:
有一个问题:它将小波系数表示为与数据形状相同的数组。
我认为将结果表示为与输入数据形状/大小相同的数组的方法是有害的。这使得整个过程变得不必要地复杂,因为你必须做出假设,或者维护一个带有索引的辅助数据结构,以便能够访问输出数组中的系数并执行逆变换(可以参考Matlab关于wavedec/waverec的文档)。
而且,尽管在理论上这个方法很好,但在实际应用中并不总是适用,因为你提到的问题:大多数情况下,输入数据的大小不是2的n次方,并且与小波滤波器卷积后的结果会比“存储空间”大,这可能导致数据丢失和不完美的重建。
为了避免这些问题,我建议使用更自然的数据结构来表示结果数据的层次,比如Python的列表、字典和元组(在适用的情况下)。