如何扩展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个边长方体.

从理论上讲,我认为应该可以确保不会出现"流血".这种小波分解和重组的算法在"数字接收器C"中讨论 - 由William Press,Saul A teukolsky,William T. Vetterling和Brian P. Flannery(第二版).虽然这种算法假定边缘处的反射而不是其他形式的边缘扩展(如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

首先,我想指出已经实现单级多维变换(源)的函数.它返回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输出相比).

基于此,将其扩展到多级案例应该相当容易.

以下是我对分解部分的看法:https://gist.github.com/934166.

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

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

将结果表示为与输入数据具有相同形状/大小的数组的方法在我看来是有害的.这让整个事情过于复杂,以理解和因为反正你必须做出的假设或保持与索引的二级数据结构能够访问系数输出数组中,并执行反变换工作(见Matlab的用于wavedec/waverec文档).

而且,即使它纸上作品伟大的,它并不总是适合,因为你刚才提到的问题,现实世界的应用:大部分的时间输入数据的大小是不是2 ^ n和小波滤波较大卷积信号的抽取结果那个"存储空间",反过来又会导致数据丢失和非完美重建.

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