PyWaveletsCWTの実装

Sep 28 2020

PyWaveletsの連続ウェーブレット変換の実装と、ここでコーディングして提供したより「基本的な」バージョンとの比較を理解しようとしています。特に：

統合されたウェーブレットは、int + diffなしで畳み込むのと同じように、畳み込み後にどのように異なりますか？
統合されたウェーブレットを、スケールを変更するのと同じ解像度でリサンプリングするにはどうすればよいですか？

ナイーブな実装：

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（関連する部分にクリップ）：

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

の比較 $f=1, 4$類洞：

回答

OverLordGoldDragon Sep 28 2020 at 20:30

まず、以下のセクションの「ナイーブな内訳」を参照してください。PyWaveletsへ：アルゴリズムはGithubで古いMATLAB実装に由来することがわかりましたが、ウェーブレットカーネルのコーディングに関する詳細は提供されていません。それで私は探検に取り掛かりました。ここにさまざまなビジュアルとコード注釈があります。主な調査結果は次のとおりです。-すべてのコード

以下の「PyWaveletsの内訳」を参照してください。それは以下を実装しています-すなわち式4：

（1）統合されたウェーブレットは、int + diffなしの畳み込みと同等で、畳み込み後にどのように異なりますか？

公式は右にあり、左から方程式を導き出します。すべてを継続的に統合する代わりに$\mathcal{R}$、積分は合計するセグメントに分割されます $k$。入力は一定であると想定されます$[k, k+1]$ （そうすべきではありませんが $)$？）、それは離散的であるのに対し、ウェーブレットは連続的に保たれます。これにより、入力を積分の外側に移動できます。視覚的に：

では、これをどのようにコーディングするのでしょうか？式4から、diffどこかに必要なことが明らかです。pywtしますnp.diff(conv)、これは基本的にnp.diff(coef); なぜナイーブの係数の差をとるように見えるのですか？

式4では、ウェーブレットを掛けk + 1たk、最大からマイナスまでの統合ウェーブレットがあります。これは、等価でウェーブレットを有する生成物へb - 1とb、それぞれ、（プラグインt=k+1とt=k、を接続して比較b=b-1し、b=binstead`）。

答えはそれからconv;にあります。ウェーブレットが信号の少なくとも半分にあるときだけになるように、すでにトリミングされていると仮定します。次に、のすべてのポイントconvは、異なるtau（またはむしろb）での入力の一部とウェーブレットの積です。だから、2つの隣接点でconvウェーブレットを持つ製品であるbとは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} $$

しかし、我々は、間にマイナスのシークt=k+1とt、以上tと下部t従ってPW内訳のステップ7を説明すると、上記の負です。

（2）スケールを変更するのと同じように、解像度を上げて統合ウェーブレットをリサンプリングするにはどうすればよいですか？

PW内訳の（5）の下の図を参照してください。pywt拡張された引数でウェーブレットを再計算するのではなく、ウェーブレットを定義するサンプルの数を介して「スケーリング」するように見えます。...しかし、視覚的に、この明確にリサンプリングすることはありません入力し、それは相対的なストレッチ、ウェーブレットを拡張します。

つまり、これは確かにスケーリングの一形態です。しかし、それはナイーブと比べてどうですか？Naiveは、固定ウェーブレット長を使用し、スケーリングされた引数を渡すことで、方程式ごとにより直感的に再計算されます。では、どちらが良いですか？反対の質問に答える方が明確です。より高いスケールとより低いスケールで何が起こるかを考えてください。

ナイーブハイ：ウェーブレットは拡張しすぎて、ゼロテールのはるか前にサンプリングフレームを終了します。
pywt high：ウェーブレットは問題ありません。問題はにありconvます。ウェーブレットの長さは1024に固定されているため、入力がこれより短い場合、スケールの大きいウェーブレットは信号を完全に乗算することはできません。視差が大きいほど、ウェーブレットは信号によって「ナイーブハイ」と同様に「見られる」ようになります。これは、垂直方向のシフトが異なる質問のヒートマップで確認できます。
- ウェーブレットの最大スケールも無制限です。len(j)スケールに比例して大きくなります（ただし、正規化の問題があります。ここを参照してください）。
ナイーブ低い：すべてが良い。
pywt lower：ウェーブレットの解像度が小さい。
固定ウィンドウ長と可変ウィンドウ長：スケールが大きいほど、入力に対するウェーブレットの拡張が大きくなるか、ゼロ以外になる必要があります。Naiveがより高いスケールで突然ゼロになることを除いて、どちらもこれを達成します。

両方のより高いスケールの問題は、最大スケールを制限することによって解決されます。両方が完全に安全なゾーンにある場合は、解像度が高いため、Naiveの実装が望ましいようです。ただし、「安全」をテールでゼロにジャンプしないと見なす場合、ナイーブのこの「安全ゾーン」は非常に小さいです。

逆に、pywtの問題は、最大ウェーブレット長を増やすことでより簡単に解決できます。これは、ジッパーのようなアーティファクト（質問のプロットで確認できます）を取り除くために提案されました。ただし、これは短い入力の大規模なパフォーマンスを悪化させますが、実際には1024より短くなる可能性はほとんどありません。

PyWaveletsの内訳：

統合前のウェーブレットは、表示されているコードblobと正確に一致します。これは、（Naiveが使用する）完全な実際のMorletの近似値です。 $\sigma > 5$Wiki。
pywtを介して実際のMorletを統合np.cumsum(psi) * stepし、異なるステップサイズを考慮します
統合されたウェーブレットはint_psi、すべての人に再利用されますscales
それぞれscaleについて、同じものint_psiが解像度を上げてリサンプリングさjれ、整数として、すべての同じ最小値と最大値（ほぼ）の間で線形に変化しますscales。

scaleナイーブと並べて（後で）増加するための結果のウェーブレット：

最初diffの畳み込み結果'd、
否定、
を掛けてsqrt(scale)、
次に、両側をトリミングして、（a）ウェーブレットが少なくとも「信号の内側」にある部分のみを含めるか、（b）その長さをに強制しlen(x)ます。

ナイーブな内訳：CWT / STFTに慣れていない場合は、ここでパート1〜3を強くお勧めします。

（1）ウェーブレットタイプが必要です。（2）ウェーブレットの長さ。（3）ウェーブレットの増分/オーバーラップ。実際のMorlet、96サンプルを使用し、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

視覚化scale=2：

問題のこれとナイーブな実装の違いは、後者の使用法はwin_inc=1、フリップされたウェーブレットで畳み込み、ウェーブレットが完全に「信号の内側」にある部分のみを含むように結果をトリミングすることと同等です。また、ここtauでは外側のループに（ "increment"）があり、スケールがありますが、2つは同等です。

最後に、タイムシフトの薄暗い方向の形状は異なります。質問の実装は、ウェーブレットが少なくとも「入力の内側」にある畳み込みの部分のみが含まれるように調整されます。これは、完全なウェーブレットが含まれる場所までさらにトリミングされます。