Wdrożenie PyWavelets CWT

Question

Staram się zrozumieć implementację ciągłej transformacji falkowej PyWavelets i jej porównanie z bardziej „podstawową” wersją, którą tutaj zakodowałem i dostarczyłem. W szczególności:

1. W jaki sposób zintegrowana falka, a następnie zróżnicowana post-konwolucja, jest równoważna splotowi bez int + diff?
2. W jaki sposób ponowne próbkowanie zintegrowanej falki przy rosnącej rozdzielczości jest równoważne zmianie jej skali?

---

Naiwna realizacja :

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 (przypięte do odpowiednich części) :

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

---

Porównanie włączone $f=1, 4$sinusoidy :

Answer 1

Najpierw zobacz „Naiwny podział” w poniższej sekcji. Na PyWavelets: algorytm został znaleziony na Github jako wywodzący się ze starej implementacji MATLAB , ale nie dostarcza żadnych szczegółów na temat kodowania jądra falkowego. Więc zacząłem odkrywać; różne wizualizacje i adnotacje do kodu tutaj . Poniżej przedstawiono kluczowe wnioski. - Cały kod

Zobacz „Podział PyWavelets” poniżej. Realizuje następujące - mianowicie równanie 4:

---

(1) W jaki sposób zintegrowana falka, a następnie zróżnicowana post-konwolucja, jest równoważna splotowi bez int + diff?

Formuła jest po prawej stronie, wyprowadzając z lewej równania; Zamiast ciągłej integracji nad wszystkim$\mathcal{R}$, całka jest dzielona na segmenty, które się sumują $k$. Zakłada się, że wejście jest stałe$[k, k+1]$ (choć nie powinno $)$?), ponieważ jest dyskretny, podczas gdy falka jest ciągła. Umożliwia to przesunięcie wejścia poza całkę; naocznie:

Jak więc to kodujemy? Z równania 4 jasno wynika, że diffgdzieś potrzebujemy . pywtrobi np.diff(conv), czyli w zasadzie np.diff(coef); dlaczego wydaje się, że bierze się pod uwagę różnicę potencjalnego współczynnika w Naiwnym?

W równaniu 4 mamy zintegrowaną falkę do k + 1, minus do k, pomnożoną przez falkę. Jest to równoważne iloczynowi z falką na poziomie b - 1i b, odpowiednio (podłącz t=k+1i t=kporównaj z podłączeniem b=b-1i b=bzamiast tego).

Odpowiedź jest już dostępna conv; załóżmy, że był już przycięty, tak że jest to tylko wtedy, gdy falka znajduje się przynajmniej w połowie wewnątrz sygnału. Wtedy każdy punkt convjest iloczynem falki z częścią wejścia, w innym tau(a raczej b). Tak więc dwa sąsiednie punkty w convto produkty z falkami w punkcie bi 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} $$

Jednak szukamy minus między t=k+1a tlub większym ti niższym t, co jest minusem powyższego, wyjaśniając w ten sposób krok 7 w Podziale PW.

---

(2) W jaki sposób ponowne próbkowanie zintegrowanej falki przy rosnącej rozdzielczości jest równoważne zmianie jej skali?

Patrz rysunek pod (5) w Podziale PW. pywtwydaje się „skalować” poprzez liczbę próbek, które definiują falkę, zamiast ponownego obliczania falki z rozszerzonym argumentem. ... choć wizualnie, to resampling wyraźnie nie rozszerzają się falki, wyciągając go w stosunku do wejścia.

Jest to więc rzeczywiście forma skalowania; ale jak to się ma do naiwności? Naiwny jest przeliczany bardziej intuicyjnie na równania, używając stałej długości falki i przekazując skalowany argument. Więc co jest lepsze? Odpowiedź na przeciwne pytanie jest jaśniejsza; zastanów się, co dzieje się w wyższej i niższej skali:

1. Naiwny wyższy : falka rozszerza się zbyt mocno, opuszczając ramkę próbkowania na długo przed zerowym ogonem.
2. pywt wyżej : falka jest w porządku; problem jest z conv. Długość falki jest ustalona na 1024, więc jeśli wejście jest krótsze, wówczas fale o większej skali nigdy nie mogą w pełni zwielokrotnić sygnału. Im większa dysproporcja, tym bardziej falka jest „widziana” podobnie do „Naiwnie wyżej” przez sygnał; widać to na mapach popularności pytania różniących się przesunięciami w pionie.
   • Maksymalna skala falki jest również nieograniczona; len(j)rośnie proporcjonalnie do skali (ale istnieje problem z normalizacją, patrz tutaj ).
3. Naiwny niższy : wszystko w porządku.
4. pywt niższy : falka ma małą rozdzielczość.
5. Stała a zmienna długość okna : wyższe skale muszą odpowiadać większemu rozszerzeniu lub niezerowej wartości falki w stosunku do wejścia . Oba osiągają to, z wyjątkiem tego, że Naive przerywa gwałtownie do zera w wyższych skalach.

W obu przypadkach problemy wyższej skali są rozwiązywane przez ograniczenie maksymalnej skali; jeśli oba są w swoich całkowicie bezpiecznych strefach, implementacja Naive wydaje się lepsza ze względu na wyższą rozdzielczość. Jednak ta „bezpieczna strefa” dla Naiwnych jest zaporowo mała, jeśli uważamy, że „bezpieczna” nie jest przeskakiwanie do zera na ogonie.

Wręcz przeciwnie, pywtproblem można łatwiej rozwiązać poprzez zwiększenie maksymalnej długości falki, co sugerowano, aby pozbyć się artefaktów przypominających zamki błyskawiczne (które można zobaczyć na wykresach pytania). To jednak pogorszyłoby wydajność w wyższej skali w przypadku krótkich danych wejściowych - ale w praktyce jest mało prawdopodobne, aby był krótszy niż 1024.

---

Podział PyWavelets :

1. Wavelet, przed integracją, pasuje dokładnie do pokazanego obiektu blob kodu, który jest przybliżeniem pełnego rzeczywistego Morleta (używanego przez Naive) przy założeniu $\sigma > 5$na Wiki .
2. pywtintegruje prawdziwy Morlet poprzez np.cumsum(psi) * step, uwzględniając różnicowy rozmiar kroku
3. Zintegrowana falka int_psijest ponownie wykorzystywana dla wszystkichscales
4. Dla każdego z nich scaleto samo int_psijest ponownie próbkowane z rosnącą rozdzielczością, przy czym jjako liczba całkowita zmienia się liniowo między tym samym minimum i maksimum (prawie) dla wszystkich scales:

1. Powstała falka do zwiększania scale, obok siebie z naiwnym (na później):

1. Wynik splotu pierwszy diff'd,
2. zanegowane,
3. pomnożone przez sqrt(scale),
4. następnie przycięte po obu stronach, aby (a) uwzględnić tylko części, w których falka znajduje się co najmniej w połowie „wewnątrz sygnału”, lub (b) wymusić jej długość len(x).

---

Naiwny podział : jeśli nie jesteś zaznajomiony z CWT / STFT, zdecydowanie polecam części 1-3 tutaj .

Potrzebujemy (1) typu falkowego; (2) długość falki; (3) przyrost / nakładanie się fali. Użyjemy prawdziwego Morleta , 96 próbek i przyrostu o 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

Wizualizowane dla scale=2:

Różnica między tą i naiwną implementacją w tym pytaniu polega na tym, że ostatnie zastosowania są win_inc=1równoważne splataniu z odwróconą falką i przycinaniu wyniku tak, aby obejmował tylko części, w których falka znajdowała się w pełni „wewnątrz sygnału”. Również tutaj mamy tau(„przyrost”) w zewnętrznej pętli, jest skala, ale te dwa są równoważne.

Wreszcie kształty wzdłuż przyciemnienia z przesunięciem czasu różnią się; implementacja pytania jest przycinana w taki sposób, że uwzględniane są tylko części splotu, w których falka jest przynajmniej w połowie „wewnątrz wejścia”; ten przycina dalej do miejsca, w którym jest uwzględniona pełna falka.

---

NIEROZWOLONE :

1. Dlaczego * sqrt(scale)(krok 8 pywt)? Eq 4 pokazuje 1 / sqrt(scale). Uważałem * scaleza rządów łańcucha, z wyjątkiem mamy zintegrowany falkowej numerycznie ... warto mieć leczenie za skalę; mamy uciec z nie zintegrowanie falki per-skali, ponieważ skala per-falkowej jest ... samo resampled falki? Rozwiązany .