PyWavelets CWT-Implementierung

Question

Ich versuche zu verstehen, wie PyWavelets die kontinuierliche Wavelet-Transformation implementiert und wie sie mit der grundlegenderen Version verglichen wird, die ich hier codiert und bereitgestellt habe. Speziell:

1. Wie ist integriertes Wavelet, das nach der Faltung differenziert wird, gleichbedeutend mit Faltung ohne int + diff?
2. Wie entspricht ein Resampling des integrierten Wavelets mit zunehmender Auflösung einer Änderung seines Maßstabs?

---

Naive Umsetzung :

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 (an relevanten Teilen befestigt) :

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

---

Vergleich auf $f=1, 4$Sinuskurven :

Answer 1

Siehe zuerst "Naive Aufschlüsselung" in einem der folgenden Abschnitte. Auf PyWavelets: Der auf Github gefundene Algorithmus stammt aus einer alten MATLAB- Implementierung , enthält jedoch keine Details zur Codierung des Wavelet-Kernels. Also machte ich mich an die Erkundung; verschiedene Grafiken und Code-Anmerkungen hier . Wichtige Erkenntnisse folgen. - Alle Code

Siehe "PyWavelets-Aufschlüsselung" weiter unten. Es implementiert Folgendes - nämlich Gleichung 4:

---

(1) Wie ist integriertes Wavelet, das nach der Faltung differenziert wird, gleichbedeutend mit Faltung ohne int + diff?

Die Formulierung ist rechts und leitet von links Gleichungen ab; Anstatt sich kontinuierlich zu integrieren$\mathcal{R}$wird das Integral in überlappende Segmente aufgeteilt $k$. Die Eingabe wird als konstant angenommen$[k, k+1]$ (sollte es aber nicht sein $)$?), da es diskret ist, während das Wavelet kontinuierlich gehalten wird. Dies ermöglicht es, den Eingang außerhalb des Integrals zu verschieben. visuell:

Wie codieren wir das? Aus Gleichung 4 geht hervor, dass wir diffirgendwo einen brauchen . pywttut np.diff(conv), was im Grunde ist np.diff(coef); Warum scheint es einen Unterschied zwischen dem möglichen Koeffizienten in Naive zu geben?

In Gleichung 4 haben wir das integrierte Wavelet bis zu k + 1minus minus kmit dem Wavelet multipliziert. Dies ist äquivalent zu Produkt mit der Wavelet bei b - 1und bjeweils (Steckern in t=k+1und t=k, zu vergleichen , mit dem Einstecken b=b-1und b=binstead`).

Die Antwort ist dann in in conv; Angenommen, es wurde bereits zugeschnitten, sodass sich das Wavelet nur dann mindestens in der Mitte des Signals befindet. Dann ist jeder Punkt in convdas Produkt des Wavelets mit einem Teil des Eingangs an einem anderen tau(oder besser gesagt b). Zwei benachbarte Punkte in convsind also Produkte mit Wavelets bei bund 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} $$

Wir suchen jedoch ein Minus zwischen t=k+1und toder ein Größeres tund Niedrigeres t, was das Negative von oben ist, und erklären so Schritt 7 in PW Breakdown.

---

(2) Wie entspricht ein Resampling des integrierten Wavelets mit zunehmender Auflösung einer Änderung seines Maßstabs?

Siehe Abbildung unter (5) in PW Breakdown. pywtscheint über die Anzahl der Abtastwerte, die das Wavelet definieren, zu "skalieren", anstatt das Wavelet mit einem erweiterten Argument neu zu berechnen. Obwohl ..., visuell, das Resampling deutlich macht die Wavelet aufzuweiten, ist es in Bezug auf den Eingang Stretching.

So it's indeed a form of scaling; but how does it compare to naive? Naive is recomputed more intuitively per equations, using fixed wavelet length, and passing a scaled argument. Then, which is better? The opposite question is clearer to answer; consider what happens at higher and lower scales:

1. Naive higher: wavelet dilates too much, exiting the sampling frame far before its zero tails.
2. pywt higher: wavelet is just fine; the problem's with conv. Wavelet length is fixed at 1024, so if the input is any shorter, then higher scale wavelets can never fully multiply the signal. The greater the disparity, the more the wavelet is "seen" similar to "Naive higher" by the signal; this can be seen in the question's heatmaps differing by vertical shifts.
   • Maximum scale of the wavelet is also unlimited; len(j) grows proportionally with scale (but there's a normalization problem, see here).
3. Naive lower: all's good.
4. pywt lower: wavelet has small resolution.
5. Fixed vs variable window length: higher scales must correspond to greater dilation, or being-nonzero, of the wavelet relative to the input. Both accomplish this, except Naive discontinues to zero abruptly at higher scales.

The higher scale problems in both are solved by limiting the maximum scale; if both are in their completely safe zones, the Naive implementation seems preferable due to higher resolution. However, this "safe zone" for Naive is prohibitively small, if we are to consider "safe" as not jumping to zero at tails.

To contrary, pywt's problem is more easily solved by increasing maximum wavelet length, which was suggested to rid of the zipper-like artifacts (which can be seen in the question's plots). That'd, however, worsen higher-scale performance for short inputs - but in practice it's unlikely to be shorter than 1024.

---

PyWavelets Breakdown:

1. Wavelet, prior to integration, matches exactly with the shown code blob, which is an approximation of the complete real Morlet (used by Naive) assuming $\sigma > 5$ in the Wiki.
2. pywt integrates real Morlet via np.cumsum(psi) * step, accounting for the differential step size
3. The integrated wavelet, int_psi, is reused for all scales
4. For each scale, the same int_psi is resampled at increasing resolution, with j as an integer varying linearly between same min and max (almost) for all scales:

1. The resulting wavelet for increasing scale, side-by-side with naive (for later):

1. The convolution result first diff'd,
2. negated,
3. multiplied by sqrt(scale),
4. then trimmed on both sides to either (a) only include parts where the wavelet is at least halfway "inside the signal", or (b) force its length to len(x).

---

Naive Breakdown: if unfamiliar with CWT/STFT, I strongly recommend parts 1-3 here.

We need (1) wavelet type; (2) wavelet length; (3) wavelet increment/overlap. We'll use real Morlet, 96 samples, and increment by 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

Visualized for scale=2:

The difference between this and naive implementation in the question is, latter uses win_inc=1, equivalent to convolving with a flipped wavelet and trimming the result to only include parts where the wavelet was fully "inside the signal". Also, here we have tau ("increment") in the outer loop, there scale, but the two are equivalent.

Lastly, shapes along timeshift dim differ; question's implementation trims so that only parts of convolution where wavelet is at least halfway "inside the input" are included; this one trims further to where the full wavelet is included.

---

UNSOLVED:

1. Why * sqrt(scale) (step 8 pywt)? Eq 4 shows 1 / sqrt(scale). I considered * scale per chain rule, except we've integrated the wavelet numerically... Makes sense to have per-scale treatment; do we get away with not integrating a per-scale wavelet because the per-scale wavelet is... the same resampled wavelet? Resolved.