テーマ： Excel によるデジタルフィルタ 観測波形をフーリエ変換してフーリエ・スペクトルを求めたのち，成分の一部を減衰さ せたスペクトルを逆フーリエ変換すると，元の観測波形にフィルタを掛けた波形が得られ る．これをデジタルフィルタという．Excel にはデータ分析ツールの一つとして，フーリエ 変換および逆フーリエ変換の機能が装備されている．この機能を利用すると手軽にデジタ ルフィルタを実現することができるが，フーリエ変換および逆変換のたびに，手動で分析 ツールの呼び出しが必要なため操作が煩雑である．自動的にデジタルフィルタリングを行 うには，分析ツールの呼び出しを VBA で自動化するか，フーリエ変換および逆変換そのも のを VBA で組む必要がある．本資料では，Excel のデータ分析ツールを手動で使用する方 法，データ分析ツールの操作を VBA で自動化する方法，VBA で FFT（高速フーリエ変換） プログラムを記述する方法について述べる． １．デジタルフィルタの概要 1.1 フーリエ変換 時刻 t における観測波形をx(t)とすると，x(t)のフーリエ変換X
は
dt t j t t x dt e t x X j t
sin cos (3) このX
は複素数で表され，実数部XR
を，虚数部をXI
と表すと
XR
jXI
X 1.2 逆フーリエ変換 フーリエ・スペクトルX
の逆変換は
dt e x X t x d e x X t x t j t j 2 1 2 1 (3) で表される． 1.3 フィルタ関数 観測波形x(t)のフーリエ変換X
に対して，フィルタ関数W
との積X
W
を計算 し，X
W
に対してフーリエ逆変換を行うと，信号成分と雑音成分を分離することが できる．フィルタ関数は，フィルタリングの種類すなわちローパス（低域通過），バンドパ ス（帯域通過），ハイパス（高域通過）によって，図 1 に示すように定義する．1.4 計算例 （１）観測波形 フーリエ変換を行うには，2 の倍数（2, 4, 8, 16, … , 2n）個のデータが必要となる．デー タ数が不足する場合は，0 を追加して 2 の倍数個にしなければならない．ここでは，図 2 に示す 1024 個のデータからなる観測波形データがあるものとする． 図 2 観測波形 （２）フーリエ・スペクトル 図 1 の観測波形にフーリエ変換を行うと図 3 に示すフーリエ・スペクトルが得られる． （ただし，N=0 のスペクトルは表示を省略）このとき，=0 の点を除き，X
は N/2 を中 心に左右対称となる． 8000 8500 9000 9500 10000 10500 11000 11500 12000 12500 1 46 91 13 6 18 1 22 6 27 1 31 6 36 1 40 6 45 1 49 6 54 1 58 6 63 1 67 6 72 1 76 6 81 1 85 6 90 1 94 6 99 1 観測値 N観測波形データ
, x(t)
1 1 1 W() W() W() N3 N4 N1 N2 N3 N4 N1 N2 低域通過フィルタ関数 帯域通過フィルタ関数 高域通過フィルタ関数 図 1 フィルタ関数図 3 フーリエ・スペクトル （３）フィルタ関数 低域通過フィルタ関数を図 4 のように定義する．このとき，=0 の点を除き，W
は N/2 を中心に左右対称となる．この例では，
1
N0,1,2,,20 and N1004,1005,,1024
W
0
N21,22,,1003
W
である． 図 4 フィルタ関数 （４）フィルタ通過後のスペクトル フィルタ通過後のスペクトルX
W
は，図 5 のとおりである． 0 100000 200000 300000 400000 500000 600000 700000 1 48 95 14 2 18 9 23 6 28 3 33 0 37 7 42 4 47 1 51 8 56 5 61 2 65 9 70 6 75 3 80 0 84 7 89 4 94 1 98 8 フ ー リ エ ・ス ペ クト ル , X ( ) Nフーリエ・スペクトル
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1 44 87 13 0 17 3 21 6 25 9 30 2 34 5 38 8 43 1 47 4 51 7 56 0 60 3 64 6 68 9 73 2 77 5 81 8 86 1 90 4 94 7 99 0 W ( ) Nフィルタ
,
W
(
)
図 5 フィルタ関数 （５）フィルタ通過後の波形 図 5 のスペクトルをフーリエ逆変換することにより，図 6 に示すように低域通過後の波 形が得られ，平滑化が行われたことが分かる． 図 6 フィルタ通過後の波形 ２．Excel を用いた FFT によるデジタルフィルタリング 2.1 Excel のデータ分析ツールを手動で使用する方法 0 100000 200000 300000 400000 500000 600000 700000 1 48 95 14 2 18 9 23 6 28 3 33 0 37 7 42 4 47 1 51 8 56 5 61 2 65 9 70 6 75 3 80 0 84 7 89 4 94 1 98 8 フ ー リ エ・ ス ペク ト ル N
フィルター通過後のスペクトル
8000 8500 9000 9500 10000 10500 11000 11500 12000 1 46 91 13 6 18 1 22 6 27 1 31 6 36 1 40 6 45 1 49 6 54 1 58 6 63 1 67 6 72 1 76 6 81 1 85 6 90 1 94 6 99 1 軸ラベル Nフィルタ通過後の波形
小西研究室フリーウェアーのページ http://www.sit.ac.jp/user/konishi/JPN/Freeware/Freeware.html より，「Excel によるデジタルフィルタ Ver. 1.0」をダウンロードし，「デジタルフィルタ(フ ーリエ解析機能利用+手動操作).xslm」という名の Excel シートを開く． 注意：ダウンロードしたシート上の観測波形のデータを単に更新しただけでは，自動計算 されません．必ず下記の要領で操作が必要になります． 操作方法①：準備 リボンのファイルタブから，オプション→アドイン→分析ツールをチェック→設定ボタン をクリックする． 「アドイン」ダイアログボックスで，「分析ツール」および「分析ツール-VBA」をチェッ クし，OK をクリックする．
操作方法②：データの更新と削除 C 列の観測波形のデータを更新する． D 列と K 列のデータを削除する． 操作方法③：データ範囲の指定 リボンのデータタブ→データ分析（下図）→フーリエ解析を選択→OK をクリック 図 の ダ イ ア ロ グ ボ ッ ク ス が 表 示 さ れ る の で ， 入 力 範 囲 を $C$4:$C$1027 ， 出 力 先 を $D$4:$D$1027 に指定し，OK をクリック 注意：D 列の値を削除しなかった場合は，次のメッセージが表示されるので，OK をクリ
ックする．D 列が空欄の場合はこのメッセージは表示されない． 操作方法④：実部，虚部，フーリエ・スペクトルの計算 E4 セルの数式を =IMREAL(D4)，F4 セルの数式を=IMAGINARY(D4)，G4 セルの数式を =IMABS(D4)とし，これらの数式を E 列 F 列 G 列の 5～1027 行に貼り付ける 操作方法⑤：フィルタ関数の指定 フィルタの値を H 列に入力する 操作方法⑥：フィルタ通過後の複素数計算 I4 セルの数式を=COMPLEX(H4*E4,H4*F4)とし，I 列の 5～1027 行に貼り付ける 操作方法⑦：フィルタ通過後のスペクトルの計算 J4 セルの数式を=IMABS(I4)とし，これらの数式を J 列の 5～1027 行に貼り付ける 操作方法⑧：逆フーリエ変換を行う リボンのデータタブ→データ分析→フーリエ解析をクリック 図 の ダ イ ア ロ グ ボ ッ ク ス が 表 示 さ れ る の で ， 入 力 範 囲 を $I$4:$I$1027 ， 出 力 先 を $K$4:$K$1027 に指定，逆変換をチェックし，OK をクリック K 列の値を削除しなかった場合は，上書き確認のメッセージが表示されるので，OK をク リックする．K 列が空欄の場合はこのメッセージは表示されない． 注意：K 列の値はテキスト形式で表示される．値が表示された直後は変換範囲がすでに選 択されているので，直ちに下図のように数値に変化する処理を行う．
以上の操作で，上述の図 6 が更新される． 2.2 データ分析ツールの操作を VBA で自動化する方法 小西研究室フリーウェアーのページ， http://www.sit.ac.jp/user/konishi/JPN/Freeware/Freeware.html より，「Excel によるデジタルフィルタ Ver. 1.0」をダウンロードし，「デジタルフィルタ(フ ーリエ解析機能利用+手動操作).xslm」という名の Excel シートを開く． 注意：ダウンロードしたシート上の観測波形のデータを単に更新しただけでは，自動計算 されません．必ず下記の要領で操作必要になります． 操作方法①：準備 Excel は直前の設定に従うため，ファイルをダウンロードしただけでは，VBA を使用でき るとは限らない．リボンに開発タブが表示されていない場合は VBA を使用するため，以下 の設定を行うこと． リボンのファイルタブから，リボンのユーザ設定を選ぶ．開発にチェックを入れ，OK をクリックしてダイアログを閉じる．これにより，メニューに「開発」が追加される．
次に，Excel の分析ツールを使用するため，2.1 の操作方法①：準備に記述した設定を行 う． 操作方法②：データの更新と削除 C 列の観測波形のデータを更新する． D 列と K 列のデータを削除する． 操作方法③：フーリエ変換 「ボタン 1」をクリックすると，フィルタ通過後の値（K 列）が自動計算される． 参考：「ボタン 1」には，以下のプログラムを記述する． Sub ボタン 1_Click() Dim i0, n As Integer Dim i0str, i1str As String i0 = Cells(1, 2)
n = Cells(2, 2)
i0str = Right$(Str(i0), Len(Str(i0) - 1))
i1str = Right$(Str(i0 + n - 1), Len(Str(i0 + n - 1) - 1))
Application.Run "ATPVBAEN.XLAM!Fourier", ActiveSheet.Range("$C$" + i0str + ":$C$" + i1str), ActiveSheet.Range("$D$" + i0str + ":$D$" + i1str), False, False
For i = i0 To i0 + n - 1
Cells(i, 5) = WorksheetFunction.ImReal(Cells(i, 4)) Cells(i, 6) = WorksheetFunction.Imaginary(Cells(i, 4)) Cells(i, 7) = WorksheetFunction.ImAbs(Cells(i, 4))
If Cells(i, 8) = 0 Then Cells(i, 9) = 0 Else
Cells(i, 9) = WorksheetFunction.Complex(Cells(i, 8) * Cells(i, 5), Cells(i, 8) * Cells(i, 6))
End If Next i
Application.Run "ATPVBAEN.XLAM!Fourier", ActiveSheet.Range("$I$" + i0str + ":$I$" + i1str), ActiveSheet.Range("$K$" + i0str + ":$K$" + i1str), True, False
For i = i0 To i0 + n - 1 Cells(i, 11) = Val(Cells(i, 11)) Next i End Sub 2.3 VBA で FFT（高速フーリエ変換）プログラムを記述する方法 小西研究室フリーウェアーのページ， http://www.sit.ac.jp/user/konishi/JPN/Freeware/Freeware.html より，「Excel によるデジタルフィルタ Ver. 1.0」をダウンロードし，「デジタルフィルタ(フ ーリエ解析機能利用+手動操作).xslm」という名の Excel シートを開く． 注意：ダウンロードしたシート上の観測波形のデータを単に更新しただけでは，自動計算 されません．必ず下記の要領で操作必要になります． 操作方法 ① Excel の VBA を使用するため，2.2 の操作方法①：準備に記述した設定を行う．ただし， VBA で FFT プログラムを記述するため，分析ツールの設定は不要． ② 観測波形のデータ（C 列）を更新する． ③ フィルタの値（G 列）を設定する． ④ ボタン 1 をクリックする． 上記操作により，フィルタ通過後の値（H 列）が自動計算される． 謝意 FFT（高速逆フーリエ変換）のプログラムは，露本伊佐男氏のホームページに記載され たものを使用しました．URL は次のとおりです． http://tsuyu.cocolog-nifty.com/blog/2007/03/publi.html 逆フーリエ変換のプログラムは，上記 FFT プログラムを逆 FFT 用に修正したものです． この場を借りてお礼申し上げます． 参考：「ボタン 1」には，以下のプログラムを記述する．
Sub ボタン 1_Click() Dim g, h, i, j, k, l, m, o, p, q As Integer i = 0 j = 0 k = 0 l = 0 p = 0 h = 0 g = 0 q = 0 'Fourier 変換 'データ数 n を指定（2 の倍数） Dim i0, n As Integer
Dim i0str, i1str As String i0 = Cells(1, 2)
n = Cells(2, 2)
i0str = Right$(Str(i0), Len(Str(i0) - 1))
i1str = Right$(Str(i0 + n - 1), Len(Str(i0 + n - 1) - 1)) 'n = 1024
m = Log(n) / Log(2)
'xr,xi はデータ数以上、s,c はデータ数の半分以上を指定しておく Dim W(1024), xr(1024), xi(1024), xd, s(512), c(512), a, b As Double 'データの読み込み For i = 1 To n xr(i - 1) = Cells(i + i0 - 1, 3) xi(i - 1) = 0 Next i 'FFT の計算 a = 0 b = 3.14159265359 * 2 / n For i = 0 To n / 2 s(i) = Sin(a) c(i) = Cos(a) a = a + b Next i l = n h = 1 For g = 1 To m l = l / 2 k = 0 For q = 1 To h p = 0 For i = k To l + k - 1
j = i + l
a = xr(i) - xr(j) b = xi(i) - xi(j) xr(i) = xr(i) + xr(j) xi(i) = xi(i) + xi(j) If p = 0 Then xr(j) = a xi(j) = b Else xr(j) = a * c(p) + b * s(p) xi(j) = b * c(p) - a * s(p) End If p = p + h Next i k = k + l + l Next q h = h + h Next g j = n / 2 For i = 1 To n - 1 k = n If j < i Then xd = xr(i) xr(i) = xr(j) xr(j) = xd xd = xi(i) xi(i) = xi(j) xi(j) = xd End If k = k / 2 Do While j >= k j = j - k k = k / 2 Loop j = j + k Next i 'データの出力 For i = 1 To n Cells(i + i0 - 1, 4) = xr(i - 1) Cells(i + i0 - 1, 5) = xi(i - 1)
Cells(i + i0 - 1, 6) = Sqr((xr(i - 1) ^ 2 + xi(i - 1) ^ 2)) Next i
'逆 Fourier 変換 'データの読み込み
For i = 1 To n
W(i - 1) = Cells(i + i0 - 1, 7) 'フィルタ関数 xr(i - 1) = xr(i - 1) * W(i - 1)
xi(i - 1) = xi(i - 1) * W(i - 1) Next i 'FFT の計算 l = n h = 1 For g = 1 To m l = l / 2 k = 0 For q = 1 To h p = 0 For i = k To l + k - 1 j = i + l a = xr(i) - xr(j) b = xi(i) - xi(j) xr(i) = xr(i) + xr(j) xi(i) = xi(i) + xi(j) If p = 0 Then xr(j) = a xi(j) = b Else xr(j) = a * c(p) + b * s(p) xi(j) = b * c(p) - a * s(p) End If p = p + h Next i k = k + l + l Next q h = h + h Next g j = n / 2 For i = 1 To n - 1 k = n If j < i Then xd = xr(i) xr(i) = xr(j) xr(j) = xd xd = xi(i) xi(i) = xi(j) xi(j) = xd End If k = k / 2 Do While j >= k
j = j - k k = k / 2 Loop j = j + k Next i Cells(i0, 8) = xr(0) / n 'データの出力 For i = 2 To n Cells(i0 + n - i + 1, 8) = xr(i - 1) / n Next i End Sub http://www.sit.ac.jp/user/konishi/JPN/Tech_inform/Pdf/DigitalFilterbyExcel.pdf
Copyright ⓒ 2013-2014 小西克享, All Rights Reserved. 個人的な学習の目的以外での使用，転載，配布等はできません．
