Implémentation PyWavelets CWT

Question

Je cherche à comprendre l'implémentation de PyWavelets de la transformation en ondelettes continue, et comment elle se compare à la version plus `` basique '' que j'ai codée et fournie ici. En particulier:

1. Comment l'ondelette intégrée, puis la post-convolution différenciée, est-elle équivalente à la convolution sans int + diff?
2. En quoi le rééchantillonnage de l'ondelette intégrée à une résolution croissante équivaut-il à changer son échelle?

---

Implémentation naïve :

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 (clipsés aux parties pertinentes) :

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

---

Comparaison sur $f=1, 4$sinusoïdes :

Answer 1

Voir d'abord "Naive Breakdown" dans une section ci-dessous. Sur PyWavelets: l'algorithme a été trouvé sur Github pour provenir d'une ancienne implémentation MATLAB , mais il ne fournit aucun détail sur le codage du noyau d'ondelettes. Alors je me suis mis à explorer; divers visuels et annotations de code ici . Les principales conclusions suivent. - Tout le code

Voir «Répartition des PyWavelets» ci-dessous. Il met en œuvre ce qui suit - à savoir Eq 4:

---

(1) Comment l'ondelette intégrée, puis la post-convolution différenciée, est-elle équivalente à la convolution sans int + diff?

La formulation est dans la droite, dérivant de la gauche, des équations; Au lieu d'intégrer en permanence sur tout$\mathcal{R}$, l'intégrale est divisée en segments ssommant sur $k$. L'entrée est supposée constante sur$[k, k+1]$ (mais ne devrait-il pas être $)$?), car elle est discrète, alors que l'ondelette est maintenue continue. Cela permet de déplacer l'entrée en dehors de l'intégrale; visuellement:

Alors, comment coder cela? Il est clair d'après l'Eq 4 que nous avons besoin d'un diffquelque part. pywtfait np.diff(conv), ce qui est fondamentalement np.diff(coef); pourquoi semble-t-il prendre la différence du coefficient potentiel dans Naive?

Dans l'Eq 4, nous avons l'ondelette intégrée up-to k + 1, moins up-to k, multipliée par l'ondelette. Ceci est équivalent au produit avec l'ondelette à b - 1et b, respectivement (branchez t=k+1et t=k, comparez à brancher b=b-1et à la b=bplace`).

La réponse est alors dedans conv; Supposons qu'il a déjà été coupé pour que ce ne soit que lorsque l'ondelette est au moins à mi-chemin dans le signal. Ensuite, chaque point convd'entrée est le produit de l'ondelette avec une partie de l'entrée, à un niveau différent tau(ou plutôt b). Donc, deux points adjacents dans convsont des produits avec des ondelettes à bet 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} $$

Cependant, nous cherchons moins entre t=k+1et t, ou un plus grand tet un plus bas t, qui est le négatif de ci-dessus, expliquant ainsi l'étape 7 dans PW Breakdown.

---

(2) En quoi le rééchantillonnage de l'ondelette intégrée à une résolution croissante équivaut-il à changer son échelle?

Voir la figure sous (5) dans PW Breakdown. pywtsemble «mettre à l'échelle» via le nombre d'échantillons qui définissent l'ondelette, plutôt que de recalculer l'ondelette avec un argument dilaté. ... bien que, visuellement, ce rééchantillonnage clairement fait dilatent l'ondelette, l' étirant par rapport à l'entrée.

C'est donc effectivement une forme de mise à l'échelle; mais comment se compare-t-il à naïf? Naive est recalculé plus intuitivement par équations, en utilisant une longueur d'ondelette fixe et en passant un argument mis à l'échelle. Alors, quel est le meilleur? La question opposée est plus claire à répondre; considérez ce qui se passe à des échelles supérieures et inférieures:

1. Naïf supérieur : l'ondelette se dilate trop, sortant de la base d'échantillonnage bien avant sa queue zéro.
2. pywt supérieur : l'ondelette est très bien; le problème est avec conv. La longueur des ondelettes est fixée à 1024, donc si l'entrée est plus courte, les ondelettes d'échelle plus élevée ne peuvent jamais multiplier complètement le signal. Plus la disparité est grande, plus l'ondelette est "vue" similaire à "Naive plus haut" par le signal; cela se voit dans les cartes thermiques de la question qui diffèrent par des décalages verticaux.
   • L'échelle maximale de l'ondelette est également illimitée; len(j)croît proportionnellement à l'échelle (mais il y a un problème de normalisation, voir ici ).
3. Naïf inférieur : tout va bien.
4. pywt lower : l'ondelette a une petite résolution.
5. Longueur de fenêtre fixe vs variable : des échelles plus élevées doivent correspondre à une plus grande dilatation, ou étant non nulle, de l'ondelette par rapport à l'entrée . Les deux accomplissent cela, sauf que Naive s'arrête brusquement à zéro à des échelles plus élevées.

Les problèmes d'échelle plus élevée dans les deux sont résolus en limitant l'échelle maximale; si les deux sont dans leurs zones complètement sûres, l'implémentation Naive semble préférable en raison d'une résolution plus élevée. Cependant, cette "zone de sécurité" pour Naive est d'une taille prohibitive, si nous devons considérer "sûre" comme ne pas sauter à zéro à la queue.

Au contraire, pywtle problème de 's est plus facilement résolu en augmentant la longueur maximale des ondelettes, ce qui a été suggéré pour se débarrasser des artefacts de type fermeture éclair (qui peuvent être vus dans les graphiques de la question). Cependant, cela aggraverait les performances à plus grande échelle pour les entrées courtes - mais en pratique, il est peu probable qu'il soit inférieur à 1024.

---

Répartition des PyWavelets :

1. L'ondelette, avant l'intégration, correspond exactement au blob de code montré, qui est une approximation du Morlet réel complet (utilisé par Naive) en supposant $\sigma > 5$dans le Wiki .
2. pywtintègre le vrai Morlet via np.cumsum(psi) * step, en tenant compte de la taille de pas différentielle
3. L'ondelette intégrée int_psi, est réutilisée pour tousscales
4. Pour chacun scale, le même int_psiest rééchantillonné à une résolution croissante, avec jcomme entier variant linéairement entre les mêmes min et max (presque) pour tous scales:

1. L'ondelette résultante pour augmenter scale, côte à côte avec naïf (pour plus tard):

1. Le résultat de la convolution en premier diff'd,
2. annulé,
3. multiplié par sqrt(scale),
4. puis rognée des deux côtés pour soit (a) inclure uniquement les parties où l'ondelette est au moins à mi-chemin «à l'intérieur du signal», soit (b) forcer sa longueur à len(x).

---

Panne naïve : si vous ne connaissez pas CWT / STFT, je recommande vivement les parties 1 à 3 ici .

Nous avons besoin (1) du type ondelette; (2) longueur d'ondelette; (3) incrément / chevauchement d'ondelettes. Nous utiliserons du vrai Morlet , 96 échantillons et un incrément de 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

Visualisé pour scale=2:

La différence entre ceci et l'implémentation naïve dans la question est, que cette dernière utilise win_inc=1, équivaut à convoluer avec une ondelette retournée et à rogner le résultat pour n'inclure que les parties où l'ondelette était complètement "à l'intérieur du signal". Aussi, ici nous avons tau("incrément") dans la boucle externe, il y a une échelle, mais les deux sont équivalents.

Enfin, les formes le long du décalage temporel diffèrent; l'implémentation de la question se règle de sorte que seules les parties de convolution où l'ondelette est au moins à mi-chemin "à l'intérieur de l'entrée" sont incluses; celui-ci est ajusté à l'endroit où l' ondelette complète est incluse.

---

NON RÉSOLU :

1. Pourquoi * sqrt(scale)(étape 8 pywt)? Eq 4 montre 1 / sqrt(scale). J'ai considéré * scalepar règle de chaîne, sauf que nous avons intégré l'ondelette numériquement ... Cela a du sens d'avoir un traitement par échelle; pouvons-nous nous en sortir sans intégrer une ondelette par échelle parce que l'ondelette par échelle est ... la même ondelette rééchantillonnée? Résolu .