Implementación de PyWavelets CWT

Question

Intento comprender la implementación de PyWavelets de la Transformada Wavelet Continua, y cómo se compara con la versión más "básica" que he codificado y proporcionado aquí. En particular:

1. ¿Cómo es la ondícula integrada, luego diferenciada después de la convolución, equivalente a la convolución sin int + diff?
2. ¿Cómo es el remuestreo de la ondícula integrada a una resolución creciente equivalente a cambiar su escala?

---

Implementación ingenua :

def cwt(x, win_len=None, win='morlet', real=False):
    N = len(x)
    win_len = win_len or N // 8

    scales = _scales(N, nv=32)
    coef = np.zeros((len(scales), N), dtype='complex128')
    kernel = morlet_kernel(win_len)
    wl2 = win_len // 2

    for i, scale in enumerate(scales):
        coef[i, :] = np.convolve(x, kernel(scale)[::-1])[wl2:-(wl2 - 1)]
    return coef

PyWavelets (recortado a las partes relevantes) :

def pywt_cwt(data, scales):
    out = np.empty((np.size(scales),) + data.shape)
    int_psi, x = integrate_wavelet(wavelet='morl', precision=10)
    
    for i, scale in enumerate(scales):
        step = x[1] - x[0]
        j = np.arange(scale * (x[-1] - x[0]) + 1) / (scale * step)
        j = j.astype(int)  # floor
        if j[-1] >= int_psi.size:
            j = np.extract(j < int_psi.size, j)
        int_psi_scale = int_psi[j][::-1]

        conv = np.convolve(data, int_psi_scale)
        coef = - np.sqrt(scale) * np.diff(conv, axis=-1)
        
        d = (coef.shape[-1] - data.shape[-1]) / 2.
        coef = coef[..., floor(d):-ceil(d)]
        out[i, ...] = coef
    return out

---

Comparación en $f=1, 4$sinusoides :

Answer 1

Primero, vea "Desglose ingenuo" en una sección a continuación. En PyWavelets: se encontró que el algoritmo en Github proviene de una implementación antigua de MATLAB , pero no proporciona detalles sobre la codificación del kernel de wavelet. Así que me puse a explorar; varios elementos visuales y anotaciones de código aquí . A continuación se presentan los hallazgos clave. - Todo el código

Consulte "Desglose de PyWavelets" a continuación. Está implementando lo siguiente, a saber, Eq 4:

---

(1) ¿Cómo es la ondícula integrada, luego diferenciada después de la convolución, equivalente a la convolución sin int + diff?

La formulación está a la derecha, derivada de la izquierda, ecuaciones; En lugar de integrarse continuamente sobre todos$\mathcal{R}$, la integral se divide en segmentos sumando $k$. La entrada se asume constante sobre$[k, k+1]$ (aunque no debería ser $)$?), ya que es discreta, mientras que la ondícula se mantiene continua. Esto permite mover la entrada fuera de la integral; visualmente:

Entonces, ¿cómo codificamos esto? Está claro por la ecuación 4 que necesitamos un difflugar. pywthace np.diff(conv), que es básicamente np.diff(coef); ¿Por qué parece tener una diferencia del posible coeficiente en Naive?

En la ecuación 4 tenemos la ondícula integrada hasta k + 1, menos hasta k, multiplicada por la ondícula. Esto es equivalente al producto con la ondícula en b - 1y b, respectivamente (conecte t=k+1y t=k, compare con conectar b=b-1y en su b=blugar`).

Entonces la respuesta está adentro conv; suponga que ya se recortó para que sea solo cuando la ondícula esté al menos a la mitad de la señal. Entonces, cada punto conves el producto de la ondícula con una porción de la entrada, en una diferente tau(o mejor dicho b). Entonces, dos puntos adyacentes en convson productos con wavelets en by b + 1:

$$ \begin{align} \text{conv} &= [s \cdot \psi_{\text{int}}(t - 0),\ s \cdot \psi_{\text{int}}(t - 1),\ ...] \\ \Rightarrow \text{DIFF}(\text{conv}) &= [s \cdot (\psi_{\text{int}}(t - 1) - \psi_{\text{int}}(t - 0)),\ ...] \end{align} $$

Sin embargo, buscamos menos entre t=k+1y t, o mayor ty menor t, que es el negativo de arriba, explicando así el paso 7 en PW Breakdown.

---

(2) ¿Cómo es el remuestreo de la ondícula integrada a una resolución creciente equivalente a cambiar su escala?

Consulte la figura bajo (5) en Desglose de PW. pywtparece "escalar" a través del número de muestras que definen la ondícula, en lugar de volver a calcular la ondícula con un argumento dilatado. ... sin embargo, visualmente, este nuevo muestreo claramente hace dilatar el tren de ondas, que se extiende en relación con la entrada.

Por tanto, es una forma de escalar; pero, ¿cómo se compara con lo ingenuo? Naive se vuelve a calcular de forma más intuitiva por ecuaciones, utilizando una longitud de onda fija y pasando un argumento escalado. Entonces, ¿cuál es mejor? La pregunta opuesta es más clara de responder; considere lo que sucede en escalas superiores e inferiores:

1. Ingenuo más alto : la ondícula se dilata demasiado y sale del marco de muestreo mucho antes de sus colas cero.
2. pywt superior : wavelet está bien; el problema está con conv. La longitud de la onda se fija en 1024, por lo que si la entrada es más corta, las ondas de mayor escala nunca pueden multiplicar completamente la señal. Cuanto mayor es la disparidad, más se "ve" la ondícula "similar" a "ingenuo más alto" por la señal; esto se puede ver en los mapas de calor de la pregunta que se diferencian por cambios verticales.
   • La escala máxima de la ondícula también es ilimitada; len(j)crece proporcionalmente con la escala (pero hay un problema de normalización, consulte aquí ).
3. Ingenuo inferior : todo está bien.
4. pywt lower : wavelet tiene una resolución pequeña.
5. Longitud de ventana fija vs variable : las escalas más altas deben corresponder a una mayor dilatación, o ser diferente de cero, de la ondícula en relación con la entrada . Ambos logran esto, excepto que Naive se interrumpe abruptamente a cero en escalas más altas.

Los problemas de mayor escala en ambos se resuelven limitando la escala máxima; si ambos están en sus zonas completamente seguras, la implementación Naive parece preferible debido a una mayor resolución. Sin embargo, esta "zona segura" para Naive es prohibitivamente pequeña, si vamos a considerar "segura" como no saltar a cero en la cola.

Por el contrario, pywtel problema se resuelve más fácilmente aumentando la longitud máxima de la ondícula, que se sugirió para eliminar los artefactos en forma de cremallera (que se pueden ver en los gráficos de la pregunta). Sin embargo, eso empeoraría el rendimiento a mayor escala para entradas cortas, pero en la práctica es poco probable que sea más corto que 1024.

---

Desglose de PyWavelets :

1. Wavelet, antes de la integración, coincide exactamente con el blob de código mostrado, que es una aproximación del Morlet real completo (utilizado por Naive) asumiendo $\sigma > 5$en el Wiki .
2. pywtintegra la vía Morlet real np.cumsum(psi) * step, teniendo en cuenta el tamaño de paso diferencial
3. La ondícula integrada`` int_psise reutiliza para todosscales
4. Para cada uno scale, int_psise vuelve a muestrear lo mismo con una resolución creciente, con jun número entero que varía linealmente entre el mismo mínimo y máximo (casi) para todos scales:

1. La ondícula resultante para aumentar scale, lado a lado con ingenuo (para más adelante):

1. El resultado de la convolución primero diff'd,
2. negado
3. multiplicado por sqrt(scale),
4. luego recortado en ambos lados para (a) solo incluir partes donde la ondícula está al menos a la mitad "dentro de la señal", o (b) forzar su longitud a len(x).

---

Desglose ingenuo : si no está familiarizado con CWT / STFT, recomiendo encarecidamente las partes 1-3 aquí .

Necesitamos (1) tipo de ondícula; (2) longitud de la ondícula; (3) incremento / superposición de ondículas. Usaremos Morlet real , 96 muestras e incrementaremos en 20.

win_len = 96  # "win" == window == wavelet
win_inc = 20
n_wins = (len(x) - win_len) // win_inc + 1

coef = np.zeros((n_wins, len(scales)), dtype='complex128')
kernel = morlet_kernel(win_len)

for tau in range(n_wins):
    start = tau * win_inc
    end   = start + win_len
    coef[tau, :] = _transform(x[start:end], kernel, scales)

def _transform(x, kernel, scales):
    coef = np.zeros(len(scales), dtype='complex128')
    for i, scale in enumerate(scales):
        psi = np.conj(kernel(scale))
        coef[i] = np.sum(x * psi / np.sqrt(scale))
    return coef

Visualizado para scale=2:

La diferencia entre esta y la implementación ingenua en la pregunta es que los últimos usos win_inc=1equivalen a convolucionar con una ondícula invertida y recortar el resultado para incluir solo las partes donde la ondícula estaba completamente "dentro de la señal". Además, aquí tenemos tau("incremento") en el bucle exterior, hay escala, pero los dos son equivalentes.

Por último, las formas a lo largo de la atenuación temporal difieren; la implementación de la pregunta se recorta para que solo se incluyan partes de la convolución donde la ondícula está al menos a la mitad "dentro de la entrada"; éste recorta aún más hasta donde se incluye la ondícula completa .

---

SIN RESOLVER :

1. ¿Por qué * sqrt(scale)(paso 8 pywt)? La ecuación 4 muestra 1 / sqrt(scale). Lo consideré * scalepor regla de cadena, excepto que hemos integrado la ondícula numéricamente ... Tiene sentido tener un tratamiento por escala; ¿Nos salimos con la suya al no integrar una ondícula por escala porque la ondícula por escala es ... la misma ondícula remuestreada? Resuelto .