¿Cómo extender pyWavelets para trabajar con datos N-dimensionales?

Question

Esto puede ser una pregunta para un foro diferente, de ser así, hágamelo saber. Noté que solo 14 personas siguen la etiqueta wavelet.

He aquí una forma elegante de extender la descomposición wavelet en pywt (paquete pyWavelets) a múltiples dimensiones. Esto debería ejecutarse fuera de la caja si pywt está instalado. La prueba 1 muestra la descomposición y recomposición de una matriz 3D. Todo, uno tiene que hacer es aumentar el número de dimensiones y el código funcionará en descomposición/recomposición con 4, 6 o incluso 18 dimensiones de datos.

He sustituido las funciones pywt.wavedec y pywt.waverec aquí. Además, en fn_dec, muestro cómo la nueva función wavedec funciona igual que la anterior.

Sin embargo, hay una trampa: representa los coeficientes wavelet como una matriz de la misma forma que los datos. Como consecuencia, con mi conocimiento limitado de wavelets, solo he podido usarlo para Haar wavelets. Otros como DB4, por ejemplo, coeficientes de sangrado sobre los límites de estos límites estrictos (no es un problema con la representación actual de los coeficientes como lista de matrices [CA, CD1 ... CDN]. Otra pega es que solo he trabajado esto con 2^N borde cuboids de datos

Teóricamente, creo que debería ser posible asegurarse de que el "sangrado" no se produce. Un algoritmo para este tipo de descomposición wavelet y recomposición se discute en "recetas numéricas en C" - por William Press, Saul A teukolsky, William T. Vetterling y Brian P. Flannery (Segunda edición). Aunque este algoritmo asume la reflexión en los bordes en lugar de las otras formas de extensiones de borde (como zpd), el método es lo suficientemente general como para trabajo para otras formas de extensión.

Cualquier sugerencia ión sobre cómo extender este trabajo a otras wavelets?

NOTA: Esta consulta también está publicada en http://groups.google.com/group/pywavelets

Gracias, 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

En primer lugar, me gustaría apuntar a la función que ya implementa Single-level Multi-dimensional Transform (Source). Devuelve un diccionario de matrices de coeficientes n-dimensionales. Los coeficientes se tratan mediante claves que describen el tipo de transformación (aproximación/detalles) aplicada a cada una de las dimensiones.

Por ejemplo para un caso 2D el resultado es un diccionario con aproximación y detalles coeficientes arrays:

>>> 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]]} 

Dónde aa es la matriz de coeficientes con aproximación transformada aplica a ambas dimensiones (LL) y da es los coeficientes matriz con transformación de detalles aplicada a la primera dimensión y transformación de aproximación aplicada a la segunda (HL) (compárese con dwt2 output).

En función de eso, debería ser bastante fácil extenderlo a la carcasa de varios niveles.

Aquí está mi opinión sobre la parte de descomposición: https://gist.github.com/934166.

También me gustaría hacer frente a un problema que menciona en su pregunta:

Hay una captura sin embargo: Se representa los coeficientes wavelet como una serie de la misma forma que los datos.

El enfoque de representar resultados como una matriz de la misma forma/tamaño que los datos de entrada es, en mi opinión, nocivo. Hace todo esto innecesariamente complejo de entender y trabajar porque de todos modos tiene que hacer suposiciones o mantener una estructura de datos secundaria con índices para poder acceder al coeficiente en la matriz de salida y realizar una transformación inversa (consulte la documentación de Matlab para wavedec/waverec)

Además, aunque funciona muy bien en papel, no siempre se adapta a las aplicaciones del mundo real debido a los problemas que ha mencionado: la mayoría de las veces el tamaño de datos de entrada no es 2^ny el resultado diezmado de la señal de convolución el filtro wavelet es más grande que el "espacio de almacenamiento", que a su vez puede conducir a la pérdida de datos y la reconstrucción no perfecta.

Para evitar estos problemas, recomendaría utilizar estructuras de datos más naturales para representar la jerarquía de datos de resultados, como las listas de Python, los diccionarios y las tuplas (donde estén disponibles).