前回、数学の記事でフーリエ変換を扱った。
本記事ではその応用例として、データに載ったノイズ(雑音)をフーリエ変換を利用して除去する方法を示す。
また、データのノイズを除去するExcelマクロを作成した。
方法
今回は離散データに対して適用する、高速フーリエ変換と呼ばれる手法を用いる。
高速フーリエ変換(FFT)とは、離散フーリエ変換(フーリエ変換の離散版のようなもの)を高速で計算するためのアルゴリズムである。
実用上のデータ処理でフーリエ変換を用いる場合は、ほぼFFTを利用すると思って良いだろう。
高速フーリエ変換の詳細については本記事では取り上げない。
高速フーリエ変換を利用したデータのノイズ除去の手順を以下に示す。
① データのサンプル数を2の累乗個にする。
② 高速フーリエ変換をする。
③ 高速フーリエ変換後のデータに窓関数を適用する。
④ 窓関数適用後のデータを逆高速フーリエ変換する。
①は、高速フーリエ変換の性質上必要となる作業である。
②は、Excelに標準装備されている高速フーリエ変換計算機能を利用する。
③で登場する窓関数とは「ある特定の区間で1、それ以外は0を取る関数」と本記事では定義する。
この窓関数をデータに適用する(掛け算する)と、その特定の区間のみのデータを抜き出すことが出来る。
今回は高速フーリエ変換後のデータにおいて、ノイズが占める領域を0として適用し、ノイズをカットする。
最後に、ノイズをカットした高速フーリエ変換後のデータを逆高速フーリエ変換することで、ノイズを除去したデータを復元する。
サンプル
今回使うサンプルは下記ボタンから入手可能である。
マクロ本体「FFT_comp.xlsm」を開くと、下のようなウィンドウが表示される。
使い方は、ノイズを除去したいデータをB列とC列に入力し、「FFT実行、ノイズ除去」ボタンをクリックするだけである。
ボタンクリック後に計算が実行され、D列にノイズ除去後のデータが自動で書き込まれる。
各グラフの表示内容は下記の通りである。
「1. 元データ」
ノイズ除去前の生データ
「2. FFT後」
生データを高速フーリエ変換した後のデータ。
「3. 窓関数適用後」
2のデータに窓関数を適用し、ノイズ成分を除去した後のデータ。
「4. 逆FFT後」
3のデータを逆高速フーリエ変換した後に得られる、ノイズを除去が完了したデータ。
2のグラフを見ると、両端には鋭いピーク、その他の領域には微小なノイズ成分が見られる。
このノイズ成分が存在する領域で0、鋭いピーク部分を1とした窓関数を適用すると3のグラフになる。
窓関数の領域はB6セルの「窓関数閾値」で設定可能であり、絶対値が閾値以上の領域では1となる。
例えば図にあるように閾値を10に設定すれば、データが10以上または10以下の値を取る場合は窓関数が1を取り、それ以外の値をノイズと判定して0を取る。
コード詳説
「FFT_comp_xlsm」には以下のマクロが保存されている。
Sub FFT()
Dim N As Long
Application.ScreenUpdating = False
Application.DisplayAlerts = False
AddIns("分析ツール").Installed = True
AddIns("分析ツール - VBA").Installed = True
N = Worksheets(1).Range(Cells(11, 2), Cells(11, 2).End(xlDown)).Rows.Count
Worksheets(1).Range("C6").Value = N
'前データ消去
Worksheets(1).Range("AA11", ActiveCell.SpecialCells(xlLastCell)).Select
Selection.Clear
'データ数
With Worksheets(1).Cells(11, 27)
.Value = 1
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 27), Cells(10 + N, 27)), _
Type:=xlFillSeries
End With
'元データを離散フーリエ変換
Application.Run "ATPVBAEN.XLAM!Fourier" _
, Worksheets(1).Range(Cells(11, 3), Cells(10 + N, 3)) _
, Worksheets(1).Range(Cells(11, 28), Cells(10 + N, 28))
'変換後のデータの実部抜出
With Worksheets(1).Cells(11, 29)
.Value = "=IMREAL(" & Cells(11, 28).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ")"
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 29), Cells(10 + N, 29))
End With
'変換後のデータの虚部抜出
With Worksheets(1).Cells(11, 30)
.Value = "=IMAGINARY(" & Cells(11, 28).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ")"
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 30), Cells(10 + N, 30))
End With
'窓関数生成
With Worksheets(1).Cells(11, 31)
.Value = "=IF(OR(ABS(" & Cells(11, 29).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ")>$B$6,ABS(" & Cells(11, 30).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ")>$B$6),1,0)"
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 31), Worksheets(1).Cells(10 + N, 31))
End With
'実部に窓関数適用
With Worksheets(1).Cells(11, 32)
.Value = "=" & Cells(11, 29).Address(RowAbsolute:=False, ColumnAbsolute:=False) & "*" & Cells(11, 31).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ""
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 32), Worksheets(1).Cells(10 + N, 32))
End With
'虚部に窓関数適用
With Worksheets(1).Cells(11, 33)
.Value = "=" & Cells(11, 30).Address(RowAbsolute:=False, ColumnAbsolute:=False) & "*" & Cells(11, 31).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ""
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 33), Worksheets(1).Cells(10 + N, 33))
End With
'複素数に変換し複素共役を取る
With Worksheets(1).Cells(11, 34)
.Value = "=IMCONJUGATE(COMPLEX(" & Cells(11, 32).Address(RowAbsolute:=False, ColumnAbsolute:=False) & "," & Cells(11, 33).Address(RowAbsolute:=False, ColumnAbsolute:=False) & "))"
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 34), Worksheets(1).Cells(10 + N, 34))
End With
'離散フーリエ変換
Application.Run "ATPVBAEN.XLAM!Fourier" _
, Worksheets(1).Range(Cells(11, 34), Cells(10 + N, 34)) _
, Worksheets(1).Range(Cells(11, 35), Cells(10 + N, 35))
'複素共役を取ってデータ数で除算
With Worksheets(1).Cells(11, 36)
.Value = "=IMCONJUGATE(" & Cells(11, 35).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ")/$C$6"
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 36), Worksheets(1).Cells(10 + N, 36))
End With
For i = i To N
Worksheets(1).Cells(10 + i, 4).Value = Worksheets(1).Cells(10 + i, 36).Value
Next
Worksheets(1).Range("A1").Select
Application.ScreenUpdating = True
Application.DisplayAlerts = True
End Subここでは全てのコードの解説はせず、本マクロで初めて用いたコードのみ見ていくことにする。
まず導入部分の
AddIns("分析ツール").Installed = True
AddIns("分析ツール - VBA").Installed = Trueで、マクロで高速フーリエ変換の機能が使用できるようになる。
実際に高速フーリエ変換を実行するコードが
Application.Run "ATPVBAEN.XLAM!Fourier" _
, Worksheets(1).Range(Cells(11, 3), Cells(10 + N, 3)) _
, Worksheets(1).Range(Cells(11, 28), Cells(10 + N, 28))である。
Application.Run “ATPVBAEN.XLAM!Fourier” , セル領域1 , セル領域2
⇒ セル領域1のデータを高速フーリエ変換し、変換後のデータをセル領域2に出力する。
また
With Worksheets(1).Cells(11, 29)
.Value = "=IMREAL(" & Cells(11, 28).Address(RowAbsolute:=False, ColumnAbsolute:=False) & ")"
.AutoFill Destination:=Worksheets(1).Range(Cells(11, 29), Cells(10 + N, 29))
End Withにある「Adress(…」はセルのアドレス(番地)を返すコードである。
セル.Address(RowAbsolute:=False, ColumnAbsolute:=False)
⇒ セルのアドレス(番地)を相対参照で返す。
セル.Address
⇒ セルのアドレス(番地)を絶対参照で返す。
(例)
Cell(1,1).Address(RowAbsolute:=False, ColumnAbsolute:=False)
⇒戻り値は「A1」
Cell(1,1).Address
⇒戻り値は「$A$1」
さらに「AutoFill Destination」はオートフィルを実行するコードである。
セルA.AutoFill Destination:=セルAを含むセル領域B
⇒ セルAを起点として、セル領域Bにオートフィルを実行する。
(例)
Range(“A1”).AutoFill Destination:=Range(“A1″:A10”)
⇒ セルA1を起点として、セルA10までオートフィルを実行する。
終わりに
離散フーリエ変換と高速フーリエ変換の詳細については、気が向いたら勉強して記事にしようと思う。
ぶっちゃけ個人的には、今はラプラス変換の復習欲の方が強い。
だが今後は自由に取れる時間が余計に限られてくる可能性が高いのが悲しい…
END
コメント