Implementasi PyWavelets CWT

Question

Saya berusaha untuk memahami implementasi PyWavelets dari Transformasi Wavelet Berkelanjutan, dan bagaimana perbandingannya dengan versi yang lebih 'dasar' yang telah saya kodekan dan berikan di sini. Khususnya:

1. Bagaimana wavelet terintegrasi, kemudian dibedakan pasca-konvolusi, setara dengan konvolusi tanpa int + diff?
2. Bagaimana resampling wavelet terintegrasi pada resolusi yang meningkat setara dengan mengubah skalanya?

---

Implementasi naif :

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 (dipotong ke bagian yang relevan) :

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

---

Perbandingan aktif $f=1, 4$sinusoid :

Answer 1

Pertama, lihat "Perincian Naif" di bagian bawah. Ke PyWavelets: algoritme ditemukan di Github untuk berasal dari implementasi MATLAB lama , tetapi tidak memberikan detail tentang pengkodean kernel wavelet. Jadi saya mulai menjelajahi; berbagai visual dan penjelasan kode di sini . Temuan kunci menyusul. - Semua kode

Lihat "Kerusakan PyWavelets" di bawah. Ini menerapkan yang berikut - yaitu Persamaan 4:

---

(1) Bagaimana wavelet terintegrasi, kemudian dibedakan pasca-konvolusi, setara dengan konvolusi tanpa int + diff?

Rumusnya ada di kanan, berasal dari kiri, persamaan; Alih-alih terus mengintegrasikan semuanya$\mathcal{R}$, integral dibagi menjadi beberapa segmen yang dihitung $k$. Input diasumsikan konstan$[k, k+1]$ (meskipun seharusnya tidak demikian $)$?), karena itu diskrit, sedangkan wavelet dipertahankan terus menerus. Ini memungkinkan untuk memindahkan masukan ke luar integral; secara visual:

Jadi bagaimana kita mengkodekan ini? Jelas dari Persamaan 4 bahwa kita membutuhkan suatu difftempat. pywttidak np.diff(conv), yang pada dasarnya np.diff(coef); mengapa tampaknya mengambil perbedaan dari koefisien calon dalam Naive?

Dalam Persamaan 4 kita memiliki wavelet terintegrasi up-to k + 1, minus up-to k, dikalikan dengan wavelet. Ini sama dengan produk dengan wavelet di b - 1dan b, masing-masing (colokkan t=k+1dan t=k, bandingkan dengan mencolokkan b=b-1dan sebagai b=bgantinya`).

Jawabannya kemudian masuk conv; asumsikan itu sudah dipangkas sehingga hanya ketika wavelet setidaknya setengah jalan di dalam sinyal. Kemudian, setiap titik masuk convadalah produk dari wavelet dengan sebagian input, pada perbedaan tau(atau lebih tepatnya b). Jadi dua titik yang berdekatanconv adalah produk dengan wavelet pada bdan 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} $$

Namun, kami mencari minus antara t=k+1dan t, atau lebih besar tdan lebih rendah t, yang merupakan negatif dari di atas, sehingga menjelaskan langkah 7 dalam Kerusakan PW.

---

(2) Bagaimana resampling wavelet terintegrasi pada resolusi yang meningkat setara dengan mengubah skalanya?

Lihat gambar di bawah (5) di Kerusakan PW. pywtmuncul untuk "menskalakan" melalui jumlah sampel yang menentukan wavelet, daripada menghitung ulang wavelet dengan argumen dilatasi. ... meskipun, secara visual, ini resampling jelas tidak membesar wavelet, peregangan relatif terhadap input.

Jadi ini memang bentuk penskalaan; tapi bagaimana jika dibandingkan dengan naif? Naive dihitung ulang secara lebih intuitif per persamaan, menggunakan panjang wavelet tetap, dan meneruskan argumen berskala. Lalu mana yang lebih baik? Pertanyaan sebaliknya lebih jelas untuk dijawab; pertimbangkan apa yang terjadi pada skala yang lebih tinggi dan lebih rendah:

1. Naive lebih tinggi : wavelet melebar terlalu banyak, keluar dari kerangka sampling jauh sebelum ekor nolnya.
2. pywt lebih tinggi : wavelet baik-baik saja; masalahnya dengan conv. Panjang wavelet ditetapkan pada 1024, jadi jika input lebih pendek, maka wavelet skala yang lebih tinggi tidak akan pernah bisa menggandakan sinyal sepenuhnya. Semakin besar disparitas, semakin banyak wavelet "terlihat" mirip dengan "Naive lebih tinggi" oleh sinyal; Hal ini dapat dilihat pada peta panas pertanyaan yang dibedakan oleh pergeseran vertikal.
   • Skala maksimum wavelet juga tidak terbatas; len(j)tumbuh secara proporsional dengan skala (tapi ada masalah normalisasi, lihat di sini ).
3. Lebih rendah naif : semuanya baik.
4. pywt lebih rendah : wavelet memiliki resolusi kecil.
5. Panjang jendela tetap vs variabel : skala yang lebih tinggi harus sesuai dengan pelebaran yang lebih besar, atau menjadi-bukan nol, dari wavelet relatif terhadap input . Keduanya mencapai ini, kecuali Naive berhenti tiba-tiba ke nol pada skala yang lebih tinggi.

Masalah skala yang lebih tinggi di keduanya diselesaikan dengan membatasi skala maksimum; jika keduanya berada di zona yang sepenuhnya aman, penerapan Naif tampaknya lebih disukai karena resolusi yang lebih tinggi. Namun, "zona aman" untuk Naive ini sangat kecil, jika kita menganggap "aman" sebagai tidak melompat ke titik nol.

Sebaliknya, pywtmasalah lebih mudah diselesaikan dengan menambah panjang wavelet maksimum, yang disarankan untuk menghilangkan artefak mirip ritsleting (yang dapat dilihat di plot pertanyaan). Namun, hal itu akan memperburuk kinerja berskala lebih tinggi untuk input pendek - tetapi dalam praktiknya tidak mungkin lebih pendek dari 1024.

---

Kerusakan PyWavelets :

1. Wavelet, sebelum integrasi, cocok persis dengan blob kode yang ditampilkan, yang merupakan perkiraan dari Morlet nyata lengkap (digunakan oleh Naive) dengan asumsi $\sigma > 5$di Wiki .
2. pywtmengintegrasikan Morlet nyata melalui np.cumsum(psi) * step, memperhitungkan ukuran langkah diferensial
3. Wavelet terintegrasi int_psi,, digunakan kembali untuk semuascales
4. Untuk setiap scale, sama int_psiyang resampled meningkatkan resolusi, dengan jsebagai integer bervariasi secara linear antara min sama dan max (hampir) untuk semua scales:

1. Gelombang yang dihasilkan untuk bertambah scale, berdampingan dengan naif (untuk nanti):

1. Hasil konvolusi pertama diff'd,
2. dinegasikan,
3. dikalikan dengan sqrt(scale),
4. kemudian dipangkas di kedua sisi untuk (a) hanya mencakup bagian di mana wavelet setidaknya setengah "di dalam sinyal", atau (b) memaksa panjangnya len(x).

---

Perincian Naif : jika tidak terbiasa dengan CWT / STFT, saya sangat merekomendasikan bagian 1-3 di sini .

Kami membutuhkan (1) jenis wavelet; (2) panjang gelombang; (3) peningkatan / tumpang tindih wavelet. Kami akan menggunakan Morlet asli , 96 sampel, dan kenaikan sebesar 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

Divisualisasikan untuk scale=2:

Perbedaan antara ini dan implementasi naif dalam pertanyaan adalah, penggunaan terakhir win_inc=1, setara dengan konvolusi dengan wavelet terbalik dan pemangkasan hasilnya hanya mencakup bagian di mana wavelet sepenuhnya "di dalam sinyal". Juga, di sini kita memiliki tau("increment") di loop luar, ada skala, tetapi keduanya setara.

Terakhir, bentuk redup timeshift berbeda; implementasi pertanyaan memangkas sehingga hanya bagian konvolusi di mana wavelet setidaknya setengah "di dalam input" yang disertakan; yang satu ini memotong lebih jauh ke tempat wavelet penuh disertakan.

---

TIDAK TERSELESAIKAN :

1. Mengapa * sqrt(scale)(langkah 8 pywt)? Persamaan 4 menunjukkan 1 / sqrt(scale). Saya mempertimbangkan * scaleper aturan rantai, kecuali kita telah mengintegrasikan wavelet secara numerik ... Masuk akal untuk memiliki perlakuan per skala; apakah kita lolos dengan tidak mengintegrasikan wavelet per skala karena wavelet per skala adalah ... wavelet resampled yang sama? Terselesaikan .