Implementazione di PyWavelets CWT

Question

Cerco di capire l'implementazione di PyWavelets della Trasformazione Wavelet continua e come si confronta con la versione più "base" che ho codificato e fornito qui. In particolare:

1. In che modo il wavelet integrato, quindi post-convolving differenziato, è equivalente al convolving senza int + diff?
2. In che modo il ricampionamento del wavelet integrato a una risoluzione crescente equivale a cambiarne la scala?

---

Implementazione 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 (ritagliati nelle parti pertinenti) :

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

---

Confronto su $f=1, 4$sinusoidi :

Answer 1

Per prima cosa vedi "Naive Breakdown" nella sezione sottostante. Su PyWavelets: l'algoritmo è stato trovato su Github per derivare da una vecchia implementazione MATLAB , ma non fornisce dettagli sulla codifica del kernel wavelet. Così ho iniziato a esplorare; varie immagini e annotazioni di codice qui . Seguono i principali risultati. - Tutto il codice

Vedi "Analisi PyWavelets" di seguito. Sta implementando quanto segue, vale a dire l'Eq 4:

---

(1) In che modo il wavelet integrato, quindi post-convolving differenziato, equivale a convolving senza int + diff?

La formulazione è a destra, derivante da sinistra, equazioni; Invece di integrarsi continuamente$\mathcal{R}$, l'integrale viene suddiviso in segmenti che si sommano $k$. Si assume che l'input sia costante$[k, k+1]$ (anche se non dovrebbe essere $)$?), poiché è discreto, mentre la wavelet è mantenuta continua. Ciò consente di spostare l'input al di fuori dell'integrale; visivamente:

Quindi come lo codifichiamo? È chiaro dall'Eq 4 che abbiamo bisogno di un diffposto da qualche parte. pywtfa np.diff(conv), che è fondamentalmente np.diff(coef); perché sembra prendere la differenza del potenziale coefficiente in Naive?

Nell'Eq 4 abbiamo il wavelet integrato up-to k + 1, meno up-to k, moltiplicato con il wavelet. Questo è equivalente al prodotto con il wavelet su b - 1e b, rispettivamente (plug in t=k+1and t=k, confronta con plug in b=b-1and b=binstead`).

La risposta è quindi dentro conv; supponiamo che sia già stato tagliato in modo che sia solo quando la wavelet è almeno a metà strada all'interno del segnale. Quindi, ogni punto in convè il prodotto del wavelet con una porzione dell'input, a un diverso tau(o meglio b). Quindi due punti adiacenticonv sono prodotti con wavelet in be 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} $$

Tuttavia, cerchiamo meno tra t=k+1e t, o un maggiore te un minore t, che è il negativo di sopra, spiegando così il passaggio 7 in PW Breakdown.

---

(2) In che modo il ricampionamento del wavelet integrato a una risoluzione crescente equivale a cambiarne la scala?

Vedere la figura sotto (5) in Ripartizione PW. pywtsembra "scalare" tramite il numero di campioni che definiscono il wavelet, piuttosto che ricalcolare il wavelet con un argomento dilatato. ... però, visivamente, questo ricampionamento chiaramente fa dilatare il wavelet, allungandolo relativo all'ingresso.

Quindi è davvero una forma di ridimensionamento; ma come ci si confronta con l'ingenuo? Naive viene ricalcolato in modo più intuitivo per equazioni, utilizzando una lunghezza wavelet fissa e passando un argomento in scala. Allora, qual è il migliore? È più chiaro rispondere alla domanda opposta; considera cosa succede a scale più alte e più basse:

1. Ingenuo più in alto : il wavelet si dilata troppo, uscendo dal frame di campionamento molto prima della sua coda zero.
2. pywt più alto : wavelet va bene; il problema è con conv. La lunghezza delle wavelet è fissata a 1024, quindi se l'ingresso è più corto, le wavelet di scala più alta non possono mai moltiplicare completamente il segnale. Maggiore è la disparità, più la wavelet viene "vista" simile a "Naive più alta" dal segnale; questo può essere visto nelle mappe di calore della domanda che differiscono per spostamenti verticali.
   • Anche la scala massima del wavelet è illimitata; len(j)cresce proporzionalmente con la scala (ma c'è un problema di normalizzazione, vedi qui ).
3. Ingenuo inferiore : va tutto bene.
4. pywt lower : wavelet ha una risoluzione piccola.
5. Lunghezza finestra fissa vs variabile : scale più alte devono corrispondere a una maggiore dilatazione, o essere diverso da zero, del wavelet rispetto all'input . Entrambi ottengono questo, tranne che Naive interrompe bruscamente a zero su scale più alte.

I problemi di scala più alta in entrambi vengono risolti limitando la scala massima; se entrambi sono nelle loro zone completamente sicure, l'implementazione Naive sembra preferibile a causa della maggiore risoluzione. Tuttavia, questa "zona sicura" per Naive è proibitivamente piccola, se vogliamo considerare "sicura" come non saltare a zero con la coda.

Al contrario, pywtil problema di è più facilmente risolvibile aumentando la lunghezza massima delle wavelet, che è stato suggerito per eliminare gli artefatti simili a cerniere (che possono essere visti nelle trame della domanda). Ciò, tuttavia, peggiorerebbe le prestazioni su scala più alta per input brevi, ma in pratica è improbabile che sia inferiore a 1024.

---

Ripartizione di PyWavelets :

1. Wavelet, prima dell'integrazione, corrisponde esattamente al blob di codice mostrato, che è un'approssimazione del Morlet reale completo (usato da Naive) assumendo $\sigma > 5$nel Wiki .
2. pywtintegra la via Morlet reale np.cumsum(psi) * step, tenendo conto della dimensione del passo differenziale
3. Il wavelet integrato,, int_psiviene riutilizzato per tuttiscales
4. Per ciascuno scale, lo stesso int_psiviene ricampionato a risoluzione crescente, con jun numero intero variabile linearmente tra lo stesso minimo e massimo (quasi) per tutti scales:

1. Il risultante wavelet per aumentare scale, fianco a fianco con ingenuo (per dopo):

1. Il risultato della convoluzione prima diff'd,
2. negato,
3. moltiplicato per sqrt(scale),
4. quindi ritagliata su entrambi i lati per (a) includere solo le parti in cui il wavelet è almeno a metà "all'interno del segnale", o (b) forzarne la lunghezza a len(x).

---

Naive Breakdown : se non hai familiarità con CWT / STFT, consiglio vivamente le parti 1-3 qui .

Abbiamo bisogno di (1) tipo wavelet; (2) lunghezza wavelet; (3) incremento / sovrapposizione wavelet. Useremo il vero Morlet , 96 campioni e incrementeremo di 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

Visualizzata per scale=2:

La differenza tra questa e l'implementazione ingenua nella domanda è, quest'ultimo utilizza win_inc=1, equivalente a convolgere con un wavelet capovolto e troncare il risultato per includere solo le parti in cui il wavelet era completamente "all'interno del segnale". Inoltre, qui abbiamo tau("incremento") nel ciclo esterno, là scala, ma i due sono equivalenti.

Infine, le forme lungo il timeshift fioco differiscono; l'implementazione della domanda viene tagliata in modo che vengano incluse solo le parti della convoluzione in cui wavelet è almeno a metà "all'interno dell'input"; questo si ritaglia ulteriormente fino al punto in cui è incluso l' intero wavelet.

---

NON RISOLTO :

1. Perché * sqrt(scale)(passaggio 8 pywt)? Eq 4 mostra 1 / sqrt(scale). Ho considerato * scalela regola della catena, tranne per il fatto che abbiamo integrato numericamente il wavelet ... Ha senso avere un trattamento per scala; ce la caviamo senza integrare un wavelet per scala perché il wavelet per scala è ... lo stesso wavelet ricampionato? Risolto .