如何扩展pyWavelets以处理N维数据？

Question

这可能是一个不同论坛的问题，如果是这样，请让我知道。我注意到只有14个人关注小波标签。

我在这里提供了一种优雅的方法，可以将pywt (pyWavelets包)中的小波分解扩展到多维。如果安装了pywt，这应该是开箱即用的。测试1显示了3D数组的分解和重组。所有需要做的就是增加维度的数量，这样代码就可以分解/重组4、6甚至18个维度的数据。

我在这里替换了pywt.wavedec和pywt.waverec函数。另外，在fn_dec中，我展示了新的wavedec函数是如何像旧函数一样工作的。

但有一个问题:它将小波系数表示为与数据形状相同的数组。因此，由于我对小波的了解有限，我只能将其用于Haar小波。其他如DB4例如在这个严格的边界的边缘出血系数(对于当前表示为数组列表CA，CD1的系数不是问题...CDN。另一个问题是，我只处理过2^N个边长方体的数据。

理论上，我认为应该可以确保“出血”不会发生。在William Press，Saul A teukolsky，William T.Vetterling和Brian P.Flannery(第二版)的"numerical decomposition in C“中讨论了这种小波分解和重新组合的算法。虽然该算法假设在边缘处反射，而不是其他形式的边缘扩展(如zpd)，但该方法足够通用，可以用于其他形式的扩展。

关于如何将这项工作扩展到其他小波，有什么建议吗？

注意:这个查询也发布在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

Answer 1

首先，我想向您指出已经实现了Single-level Multi-dimensional Transform (Source)的函数。它返回n维系数数组的字典。系数由描述应用于每个维度的变换类型(近似/细节)的关键字来寻址。

例如，对于2D情况，结果是具有近似系数数组和详细系数数组的字典：

>>> 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 output相比)。

在此基础上，将其扩展到多级情况应该是相当容易的。

下面是我对分解部分的看法：.

我还想谈谈你在问题中提到的一个问题：

有一个问题:它将小波系数表示为与数据形状相同的数组。

在我看来，将结果表示为与输入数据形状/大小相同的数组的方法是有害的。这使得理解和处理整个事情变得不必要地复杂，因为不管怎样，您必须做出假设或维护一个带有索引的辅助数据结构，以便能够访问输出数组中的系数并执行逆转换(请参阅Matlab的文档了解wavedec/waverec)。

此外，尽管它在纸面上效果很好，但由于您提到的问题，它并不总是适合实际应用:大多数情况下，输入数据大小不是2^n，并且使用小波滤波器卷积信号的抽取结果大于“存储空间”，这反过来可能导致数据丢失和不完美的重建。

为了避免这些问题，我建议使用更自然的数据结构来表示结果数据层次结构，比如Python的列表、字典和元组(如果有)。