つれづれなる備忘録

日々の発見をあるがままに綴る

Pythonによるデータ処理10 ~ ウェーブレット変換2

 今回はPythonのPywaveletsモジュールを用いたウェーブレット変換について紹介する。

1. Pywaveletモジュール

前回はScipyのcwt関数を用いた連続ウェーブレット変換について紹介したが

atatat.hatenablog.com

今回はPywaveletsモジュール

PyWavelets - Wavelet Transforms in Python — PyWavelets Documentation

を用いたウェーブレット変換を行う。Scipyとの差は、基底関数の種類、離散ウェーブレット変換、2Dウェーブレット変換などウェーブレット変換に関して多くの関数が利用できる点にある。

Scipyと比較するためウェーブレット変換を適用する信号は前回と同じとする。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline

t = np.linspace(-1, 1, 200, endpoint=False)
tri = np.cos(2 * np.pi * 7 * t)
pulse = signal.gausspulse(t - 0.4, fc=2)
noise = np.random.normal(0,0.5,np.size(t))

sig  =  tri+ pulse +noise 

plt.figure(figsize=(8,6))
plt.plot(t,sig)
plt.plot(t,tri-3,'r')
plt.plot(t,pulse-6,'b')
plt.plot(t,noise-8,'g')
plt.legend(['sig','tri','pulse','noise'])

"テスト信号"
テスト信号

2. 連続ウェーブレット変換

 Pywaveletsモジュールがインストールされていれば、import pywtとすればロードできる。連続ウェーブレット変換はScipyと同名のcwt関数を使用する。引数の順番がことなるが、Scipyとほぼ同じ使い勝手でpywt.cwt('sig, width, 'mexh')とすると基底関数をメキシカンハット関数(Scipyではricker)とするウェーブレット変換を実行できる。widthについてはScipyと同じで1から30のアレイで指定する。

import pywt

width = np.arange(1, 31)
cwtmatr,freq = pywt.cwt(sig, width, 'mexh')
plt.imshow(cwtmatr,aspect='auto')
plt.colorbar()

連続ウェーブレット変換の結果もScipyと同じになる。

"Pywaveletのcwt関数"
Pywaveletのcwt関数

Pywaveletsモジュールでサポートされている基底関数を以下のように調べることができる。wavelist関数の引数kindは連続ウェーブレットではcontinuous、離散ウェーブレットではdiscreteを指定する。(両方の場合はall)

wavlist = pywt.wavelist(kind='continuous')
wavlist
>['cgau1',
 'cgau2',
 'cgau3',
 'cgau4',
 'cgau5',
 'cgau6',
 'cgau7',
 'cgau8',
 'cmor',
 'fbsp',
 'gaus1',
 'gaus2',
 'gaus3',
 'gaus4',
 'gaus5',
 'gaus6',
 'gaus7',
 'gaus8',
 'mexh',
 'morl',
 'shan']

試しにリストにあるShannon関数shanを用いてウェーブレット変換しようとすると、帯域幅と中心周波数を関数名で指定するように警告が表示される。

width = np.arange(1, 31)
cwtmatr,freq = pywt.cwt(sig, width, 'shan')
plt.imshow(np.abs(cwtmatr),aspect='auto')
plt.colorbar()
>/usr/local/lib/python3.7/dist-packages/pywt/_cwt.py:117: FutureWarning: Wavelets from the family shan, without parameters specified in the name are deprecated. The name should take the form shanB-C where B and C are floats representing the bandwidth frequency and center frequency, respectively (example: shan1.5-1.0).
  wavelet = DiscreteContinuousWavelet(wavelet)

警告の例の通りにshan1.5-1.0とすると、警告なしで実行される。

width = np.arange(1, 31)
cwtmatr,freq = pywt.cwt(sig, width, 'shan1.5-1.0')
plt.imshow(np.abs(cwtmatr),aspect='auto')
plt.colorbar()

"Shannon関数によるウェーブレット変換"
Shannon関数によるウェーブレット変換

3. 離散ウェーブレット変換と逆変換

 Pywaveletsモジュールを使用すると離散ウェーブレット変換をdwt関数を用いてdwt(sig.'db1')を実行し係数を得ることができる。ここでdb1は基底関数でどのような種類があるかはpywt.wavelist(kind='discrete')で確認できる。

width = np.arange(1, 31)
cA,cD = pywt.dwt(sig,  'db1')
plt.figure(figsize=(10,5))
plt.plot(cA)
plt.plot(cD)

"離散ウェーブレット変換"
離散ウェーブレット変換

フーリエ変換と同様に逆変換により元の信号に戻すことができる。離散ウェーブレット変換により得られた係数cA,cDを用いて離散ウェーブレット逆変換をidwt関数で実行すると、信号を復元することができる。

sig_r1 = pywt.idwt(cA, None ,'db1','smooth')
sig_r2 = pywt.idwt(None, cD ,'db1','smooth')
plt.figure(figsize=(10,5))
plt.plot(sig_r1+sig_r2)
plt.plot(sig,'r--')
plt.legend(['IDWT','Original'])

"離散逆ウェーブレット変換"
離散逆ウェーブレット変換

4. まとめ

 今回はPywaveletsモジュールを用いた連続ウェーブレット変換、離散ウェーブレット変換について紹介した。

Pythonによるデータ処理9 ~ ウェーブレット変換

 今回はPythonによるウェーブレット変換について紹介する。

1. ウェーブレット変換

 ウェーブレット変換はフーリエ変換のような周波数解析の1種で基底関数としてウェーブレット関数Ψ (t)を用いて、連続ウェーブレット変換は以下のように積分変換する。

 \displaystyle X (a,b) = \frac{1}{a} \int^{\infty}_{-\infty} x(t) \Psi^{*} (\frac{t-b}{a} ) dt

基本的には信号がウェーブレット関数の拡大・縮小、平行移動によってあらわす。ウェーブレット関数は許容条件がありメキシカンハット関数などが使われる。

2. 信号生成

ウェーブレット変換の例として、以下のサンプルコードを基にテスト用の信号を作成する。

scipy.signal.cwt — SciPy v1.7.1 Manual

テスト用信号sigは、cos波tri、ガウシアンパルスpulse、白色ノイズnoiseを加算したものとする。

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
%matplotlib inline

t = np.linspace(-1, 1, 200, endpoint=False)
tri = np.cos(2 * np.pi * 7 * t)
pulse = signal.gausspulse(t - 0.4, fc=2)
noise = np.random.normal(0,0.5,np.size(t))

sig  =  tri+ pulse +noise 

plt.figure(figsize=(8,6))
plt.plot(t,sig)
plt.plot(t,tri-3,'r')
plt.plot(t,pulse-6,'b')
plt.plot(t,noise-8,'g')
plt.legend(['sig','tri','pulse','noise'])

上記のコードを実行して、テスト用の信号、cos波、ガウシアンパルス、ノイズをそれぞれ表示した。

"テスト信号"
テスト信号

3. 連続ウェーブレット変換

 Pythonのウェーブレット変換できるモジュールはいくつかあるが、今回はScipyのcwt関数を利用する。使い方としては上記で生成したテスト用信号sigとしてcwtmatr = signal.cwt(sig, signal.ricker, width)とするとウェーブレット関数としてricker (メキシカンハット関数)、スケーリングパラメータ:widthで連続ウェーブレット変換を実行する。widthはアレイで生成し、np.arange(1,31)として1から30までの値としている。

widths = np.arange(1, 31)
cwtmatr = signal.cwt(sig, signal.ricker, widths)
plt.imshow(cwtmatr,aspect='auto')
plt.colorbar()

cwt関数を実行して、imshowを用いてマップ表示する。スケーリング係数が5ぐらいまではテスト用信号に含まれるcos波成分、スケーリング係数10ぐらいはガウシアンパルス成分をよくあらわしている。

"ウェーブレット変換1"
ウェーブレット変換1

Scipyのウェーブレット基底用関数としてricker以外にmorletも使用できる。ただし複素数で出力されるので、実数に変換して表示する必要がある。

cwtmatr = signal.cwt(sig, signal.morlet, widths)
plt.imshow(np.abs(cwtmatr),aspect='auto')
plt.colorbar()

"ウェーブレット変換, morlet関数"
ウェーブレット変換, morlet関数

ウェーブレット変換に使用したRicker関数を直接プロットするとスケーリング係数と関数の形の関係がわかる。ricker(points, a)のpointsはプロット点数、aはスケーリングパラメータになり、a=1だと幅が狭、a=30だと幅が広いメキシカンハット関数になっていることが確認できる。

plt.plot(signal.ricker(100,1),'g')
plt.plot(signal.ricker(100,10),'r')
plt.plot(signal.ricker(100,30),'b')

"Ricker関数"
Ricker関数

Ricker関数と同じようにプロット点数とスケーリングパラメータを引数とするユーザ定義関数をウェーブレット変換に使用することができる。以下は指数関数をウェーブレット変換に使用したが、指数関数自体は基底関数の許容条件を満たしておらず、単にユーザ定義関数がcwt関数によるウェーブレット変換に使えるということのみを示している。   

def mexp(points,a):
  x = np.linspace(-1,1,points)
  y = np.exp(-x**2/(a/50)**2)
  return y

plt.figure(figsize=(12,4))

plt.subplot(1,2,1)
cwtmatr = signal.cwt(sig, mexp, widths)
plt.imshow(cwtmatr,aspect='auto')
plt.colorbar()
plt.subplot(1,2,2)
plt.plot(mexp(200,1),'g')
plt.plot(mexp(200,10),'r')
plt.plot(mexp(200,30),'b')

"自作関数によるウェーブレット変換"
自作関数によるウェーブレット変換

3. まとめ

 今回はScipyモジュールのcwt関数による連続ウェーブレット変換の方法を紹介した。

ブログデザイン備忘録 ~2022年のトレンドカラー

今年最初の記事は、2022年のWebデザインのトレンドカラーについて取り上げたい。

coliss.com

によると2022年の色は「最も幸せで暖かいブルー」ベリーペリ(Very Peri)と名付けられ、カラーコード:#6667ab 具体的には以下のラベンダー系の色になる。

早速Very Periを使って背景デザインを行ってみた。元のデザインは

サイトの引き立て役はコレ! おしゃれすぎる背景をコピペで実装 【 HTML/CSS 】 - デシノン

の"ふわふわ四角"でHTMLとCSSのみでアニメーション背景が作成できる。HTMLは以下の通りで、基本的にはlist構造の部分を利用して四角形が回転移動するアニメーションを作成する。

<div class="context">
        <h1>Color of the year 2022</h1>
    </div>

<div class="area" >
            <ul class="circles">
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
                    <li></li>
            </ul>
    </div >

CSSでは背景をVery Periにするためbackground:#6667ab;に変更する。四角形はlist-style:none;としてマーカーをなしにして代わりにbackgound: rgba(255, 255, 255, 0.3);とすることで背景を白の透過色に設定、サイズはwidthとheightで設定、位置はleftで設定できる。10個のリストに対してそれぞれサイズ、位置を適当に設定することで様々な大きさの四角形を表示できる。 縦の移動と回転については@keyframes animate内のtransform: translateY(-1000px) rotate(720deg);として、タイミング・速度はanimation-delayanimation-durationで調整できる。

.area{
    background:#6667ab;  
    width: 100%;
    height:100vh;   
}

.circles li{
    position: absolute;
    display: block;
    list-style: none;
    width: 20px;
    height: 20px;
    background: rgba(255, 255, 255, 0.3);
    animation: animate 25s linear infinite;
    bottom: -150px;
    
}

.circles li:nth-child(1){
    left: 25%;
    width: 80px;
    height: 80px;
    animation-delay: 0s;
}

/* 中略 */
@keyframes animate {

    0%{
        transform: translateY(0) rotate(0deg);
        opacity: 1;
        border-radius: 0;
    }

    100%{
        transform: translateY(-1000px) rotate(720deg);
        opacity: 0;
        border-radius: 50%;
    }

}

Code Penの実行例を以下の通り。

See the Pen VeryPeri by ATATAT (@atatat) on CodePen.