kaggleのkernel(今はnotebook)で「こんな方法があるんだー」と納得したノイズ除去法があったので、あれこれ調べてみたことを以下にメモがてらまとめます。
##方法の概要
ある信号や波形に対しての解析手法としてウェーブレット変換と呼ばれるものがある。
解析手法としてフーリエ変換をイメージすることが多いと思うが、フーリエ変換が波形を三角関数(さいん、こさいん)の和に展開するのに対して、ウェーブレット変換では好きな形の波形で分解することが出来る。(しかも逆変換ある)
よって、与えられた波形に一部特徴的な波(マザーウェーブレット)が確認できたら、その波を基底として分解してあげることにより、どの時間帯や規則でその特徴が現れるかを確認することが出来る。(分解したときに現れる係数をウェーブレット係数という)
したがって、ノイズの形が把握出来たらそれをマザーウェーブレットとしてウェーブレット変換をし、ノイズ部分を0に置き換えてから逆変換を行えば、ノイズが無い波形が得られる。
##手順
さて、以上の仕組みを与えた関数が以下である
import numpy as np
import pywt
def maddest(d, axis=None):
return np.mean(np.absolute(d - np.mean(d, axis)), axis)
def denoise(x, wavelet='db', level=1):
coeff = pywt.wavedec(x, wavelet, mode="per")
sigma = (1/0.6745) * maddest(coeff[-level])
uthresh = sigma * np.sqrt(2*np.log(len(x)))
coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])
return pywt.waverec(coeff, wavelet, mode='per')
一番目のdefでは平均絶対偏差の計算を与えている
偏差とは各データと平均との差であり、絶対偏差とは偏差の絶対値、そして平均絶対偏差とは絶対偏差を母数で割ったもの、つまり平均である。
例:50点、70点、90点、70点の4つのデータ(平均70点)
偏差は-20、0、20、0 であり、
絶対偏差は20、0、20、0 であり 、
平均絶対偏差は(0+20+20+0)/4=10である。
二番目のdefでは波形のノイズ除去を与えている。まず
coeff = pywt.wavedec(x, wavelet, mode='per')
では与えられた波形をウェーブレット変換している。
パラメータxは配列データ、waveletはウェーブレットの名前(今回はdb)、modeは境界条件のこと(perで波形の右端と左端が繋ぐことを意味する)
coeffはウェーブレット係数のリストとして与えられる。(但し、初めから一番目はウェーブレット係数ではない。最小の解像度の余りが格納されている。)
次に
sigma = (1/0.6745) * maddest(coeff[-level])
ではノイズの大きさを与えている。maddest(coeff[-level])でウェーブレット係数のリストの後ろから一番目の値の平均絶対偏差を表している(リストの後ろであればあるほど解像度が高い。その設定をlevelとしている)。0.6745はおまじない、、といいたいところではあるが、標準偏差と平均絶対偏差との比率を表す定数である。
次に
uthresh = sigma * np.sqrt(2*np.log(len(x)))
coeff[1:] = (pywt.threshold(i, value=uthresh, mode='hard') for i in coeff[1:])
では、まず閾値(uthresh)の計算をして、ウェーブレット係数がその閾値を下回ったらノイズとみなして0を与える処理を行っている。
最後に
return pywt.waverec(coeff, wavelet, mode='per')
で、ウェーブレット変換の逆変換で元の形に戻す。先ほどの閾値とのギャップを調べる処理を経ているから、ノイズがなくなった波形が得られるのである。
##感想
もうちょっと理論おっかけてみたいと感じた。
##参考
Discrete Wavelet Transform (DWT)
ウェーブレット変換の基本
ウェーブレットによるノイズ除去