埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-1/14 テーマ: Excel によるデジタルフィルタ 観 測 波 形 をフーリエ 変 換 してフーリエ スペクトルを 求 めたのち, 成 分 の 一 部 を 減 衰 さ せたスペクトルを 逆 フーリエ 変 換 すると, 元 の 観 測 波 形 にフィルタを 掛 けた 波 形 が 得 られ る.これをデジタルフィルタという.Excel にはデータ 分 析 ツールの 一 つとして,フーリエ 変 換 および 逆 フーリエ 変 換 の 機 能 が 装 備 されている.この 機 能 を 利 用 すると 手 軽 にデジタ ルフィルタを 実 現 することができるが,フーリエ 変 換 および 逆 変 換 のたびに, 手 動 で 分 析 ツールの 呼 び 出 しが 必 要 なため 操 作 が 煩 雑 である. 自 動 的 にデジタルフィルタリングを 行 うには, 分 析 ツールの 呼 び 出 しを VBA で 自 動 化 するか,フーリエ 変 換 および 逆 変 換 そのも のを VBA で 組 む 必 要 がある. 本 資 料 では,Excel のデータ 分 析 ツールを 手 動 で 使 用 する 方 法,データ 分 析 ツールの 操 作 を VBA で 自 動 化 する 方 法,VBA で FFT( 高 速 フーリエ 変 換 ) プログラムを 記 述 する 方 法 について 述 べる. 1.デジタルフィルタの 概 要 1.1 フーリエ 変 換 時 刻 t における 観 測 波 形 を x(t) とすると, (t) この X x t e jt dt x のフーリエ 変 換 X xt cost j sin tdt X は 複 素 数 で 表 され, 実 数 部 XR を, 虚 数 部 を XI X X R jx I は と 表 すと (3) 1.2 逆 フーリエ 変 換 フーリエ スペクトル x x t で 表 される. 1 X 2 1 2 X jt jt t X x e dt x e の 逆 変 換 は d (3) 1.3 フィルタ 関 数 観 測 波 形 x(t) のフーリエ 変 換 X に 対 して,フィルタ 関 数 W との 積 XW を 計 算 X W に 対 してフーリエ 逆 変 換 を 行 うと, 信 号 成 分 と 雑 音 成 分 を 分 離 することが し, できる.フィルタ 関 数 は,フィルタリングの 種 類 すなわちローパス( 低 域 通 過 ),バンドパ ス( 帯 域 通 過 ),ハイパス( 高 域 通 過 )によって, 図 1 に 示 すように 定 義 する.
1 46 91 136 181 226 271 316 361 406 451 496 541 586 631 676 721 766 811 856 901 946 991 観 測 値 埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-2/14 W() W() W() 1 1 1 N3 N4 N1 N2 N3 N4 N1 N2 低 域 通 過 フィルタ 関 数 帯 域 通 過 フィルタ 関 数 高 域 通 過 フィルタ 関 数 図 1 フィルタ 関 数 1.4 計 算 例 (1) 観 測 波 形 フーリエ 変 換 を 行 うには,2 の 倍 数 (2, 4, 8, 16,, 2 n ) 個 のデータが 必 要 となる.デー タ 数 が 不 足 する 場 合 は,0 を 追 加 して 2 の 倍 数 個 にしなければならない.ここでは, 図 2 に 示 す 1024 個 のデータからなる 観 測 波 形 データがあるものとする. 12500 12000 11500 11000 10500 10000 9500 9000 8500 8000 観 測 波 形 データ, x(t) N 図 2 観 測 波 形 (2)フーリエ スペクトル 図 1 の 観 測 波 形 にフーリエ 変 換 を 行 うと 図 3 に 示 すフーリエ スペクトルが 得 られる. (ただし,N=0 のスペクトルは 表 示 を 省 略 )このとき,=0 の 点 を 除 き, X は N/2 を 中 心 に 左 右 対 称 となる.
1 44 87 130 173 216 259 302 345 388 431 474 517 560 603 646 689 732 775 818 861 904 947 990 W() 1 48 95 142 189 236 283 330 377 424 471 518 565 612 659 706 753 800 847 894 941 988 フーリエ スペクトル, X() 埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-3/14 700000 フーリエ スペクトル 600000 500000 400000 300000 200000 100000 0 N 図 3 フーリエ スペクトル (3)フィルタ 関 数 低 域 通 過 フィルタ 関 数 を 図 4 のように 定 義 する.このとき,=0 の 点 を 除 き, W は N/2 を 中 心 に 左 右 対 称 となる.この 例 では, W 1 N 0,1,2,,20 and N 1004,1005,,1024 W 0 N 21,22,,1003 である. 1.2 フィルタ,W() 1.0 0.8 0.6 0.4 0.2 0.0 N 図 4 フィルタ 関 数 (4)フィルタ 通 過 後 のスペクトル フィルタ 通 過 後 のスペクトル X W は, 図 5 のとおりである.
1 46 91 136 181 226 271 316 361 406 451 496 541 586 631 676 721 766 811 856 901 946 991 軸 ラベル 1 48 95 142 189 236 283 330 377 424 471 518 565 612 659 706 753 800 847 894 941 988 フーリエ スペクトル 埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-4/14 700000 600000 フィルター 通 過 後 のスペクトル 500000 400000 300000 200000 100000 0 N 図 5 フィルタ 関 数 (5)フィルタ 通 過 後 の 波 形 図 5 のスペクトルをフーリエ 逆 変 換 することにより, 図 6 に 示 すように 低 域 通 過 後 の 波 形 が 得 られ, 平 滑 化 が 行 われたことが 分 かる. 12000 11500 11000 10500 フィルタ 通 過 後 の 波 形 10000 9500 9000 8500 8000 N 図 6 フィルタ 通 過 後 の 波 形 2.Excel を 用 いた FFT によるデジタルフィルタリング 2.1 Excel のデータ 分 析 ツールを 手 動 で 使 用 する 方 法
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-5/14 小 西 研 究 室 フリーウェアーのページ http://www.sit.ac.jp/user/konishi/jpn/freeware/freeware.html より, Excel によるデジタルフィルタ Ver. 1.0 をダウンロードし, デジタルフィルタ(フ ーリエ 解 析 機 能 利 用 + 手 動 操 作 ).xslm という 名 の Excel シートを 開 く. 注 意 :ダウンロードしたシート 上 の 観 測 波 形 のデータを 単 に 更 新 しただけでは, 自 動 計 算 されません. 必 ず 下 記 の 要 領 で 操 作 が 必 要 になります. 操 作 方 法 1: 準 備 リボンのファイルタブから,オプション アドイン 分 析 ツールをチェック 設 定 ボタン をクリックする. アドイン ダイアログボックスで, 分 析 ツール および 分 析 ツール-VBA をチェッ クし,OK をクリックする.
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-6/14 操 作 方 法 2:データの 更 新 と 削 除 C 列 の 観 測 波 形 のデータを 更 新 する. D 列 と K 列 のデータを 削 除 する. 操 作 方 法 3:データ 範 囲 の 指 定 リボンのデータタブ データ 分 析 ( 下 図 ) フーリエ 解 析 を 選 択 OK をクリック 図 の ダ イ ア ロ グ ボ ッ ク ス が 表 示 さ れ る の で, 入 力 範 囲 を $C$4:$C$1027, 出 力 先 を $D$4:$D$1027 に 指 定 し,OK をクリック 注 意 :D 列 の 値 を 削 除 しなかった 場 合 は, 次 のメッセージが 表 示 されるので,OK をクリ
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-7/14 ックする.D 列 が 空 欄 の 場 合 はこのメッセージは 表 示 されない. 操 作 方 法 4: 実 部, 虚 部,フーリエ スペクトルの 計 算 E4 セルの 数 式 を =IMREAL(D4),F4 セルの 数 式 を=IMAGINARY(D4),G4 セルの 数 式 を =IMABS(D4)とし,これらの 数 式 を E 列 F 列 G 列 の 5~1027 行 に 貼 り 付 ける 操 作 方 法 5:フィルタ 関 数 の 指 定 フィルタの 値 を H 列 に 入 力 する 操 作 方 法 6:フィルタ 通 過 後 の 複 素 数 計 算 I4 セルの 数 式 を=COMPLEX(H4*E4,H4*F4)とし,I 列 の 5~1027 行 に 貼 り 付 ける 操 作 方 法 7:フィルタ 通 過 後 のスペクトルの 計 算 J4 セルの 数 式 を=IMABS(I4)とし,これらの 数 式 を J 列 の 5~1027 行 に 貼 り 付 ける 操 作 方 法 8: 逆 フーリエ 変 換 を 行 う リボンのデータタブ データ 分 析 フーリエ 解 析 をクリック 図 の ダ イ ア ロ グ ボ ッ ク ス が 表 示 さ れ る の で, 入 力 範 囲 を $I$4:$I$1027, 出 力 先 を $K$4:$K$1027 に 指 定, 逆 変 換 をチェックし,OK をクリック K 列 の 値 を 削 除 しなかった 場 合 は, 上 書 き 確 認 のメッセージが 表 示 されるので,OK をク リックする.K 列 が 空 欄 の 場 合 はこのメッセージは 表 示 されない. 注 意 :K 列 の 値 はテキスト 形 式 で 表 示 される. 値 が 表 示 された 直 後 は 変 換 範 囲 がすでに 選 択 されているので, 直 ちに 下 図 のように 数 値 に 変 化 する 処 理 を 行 う.
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-8/14 以 上 の 操 作 で, 上 述 の 図 6 が 更 新 される. 2.2 データ 分 析 ツールの 操 作 を VBA で 自 動 化 する 方 法 小 西 研 究 室 フリーウェアーのページ, http://www.sit.ac.jp/user/konishi/jpn/freeware/freeware.html より, Excel によるデジタルフィルタ Ver. 1.0 をダウンロードし, デジタルフィルタ(フ ーリエ 解 析 機 能 利 用 + 手 動 操 作 ).xslm という 名 の Excel シートを 開 く. 注 意 :ダウンロードしたシート 上 の 観 測 波 形 のデータを 単 に 更 新 しただけでは, 自 動 計 算 されません. 必 ず 下 記 の 要 領 で 操 作 必 要 になります. 操 作 方 法 1: 準 備 Excel は 直 前 の 設 定 に 従 うため,ファイルをダウンロードしただけでは,VBA を 使 用 でき るとは 限 らない.リボンに 開 発 タブが 表 示 されていない 場 合 は VBA を 使 用 するため, 以 下 の 設 定 を 行 うこと. リボンのファイルタブから,リボンのユーザ 設 定 を 選 ぶ. 開 発 にチェックを 入 れ,OK をクリックしてダイアログを 閉 じる.これにより,メニューに 開 発 が 追 加 される.
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-9/14 う. 次 に,Excel の 分 析 ツールを 使 用 するため,2.1 の 操 作 方 法 1: 準 備 に 記 述 した 設 定 を 行 操 作 方 法 2:データの 更 新 と 削 除 C 列 の 観 測 波 形 のデータを 更 新 する. D 列 と K 列 のデータを 削 除 する. 操 作 方 法 3:フーリエ 変 換 ボタン 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))
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-10/14 6)) 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, End If 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)) End Sub 2.3 VBA で FFT( 高 速 フーリエ 変 換 )プログラムを 記 述 する 方 法 小 西 研 究 室 フリーウェアーのページ, http://www.sit.ac.jp/user/konishi/jpn/freeware/freeware.html より, Excel によるデジタルフィルタ Ver. 1.0 をダウンロードし, デジタルフィルタ(フ ーリエ 解 析 機 能 利 用 + 手 動 操 作 ).xslm という 名 の Excel シートを 開 く. 注 意 :ダウンロードしたシート 上 の 観 測 波 形 のデータを 単 に 更 新 しただけでは, 自 動 計 算 されません. 必 ず 下 記 の 要 領 で 操 作 必 要 になります. 操 作 方 法 1 Excel の VBA を 使 用 するため,2.2 の 操 作 方 法 1: 準 備 に 記 述 した 設 定 を 行 う.ただし, VBA で FFT プログラムを 記 述 するため, 分 析 ツールの 設 定 は 不 要. 2 観 測 波 形 のデータ(C 列 )を 更 新 する. 3 フィルタの 値 (G 列 )を 設 定 する. 4 ボタン 1 をクリックする. 上 記 操 作 により,フィルタ 通 過 後 の 値 (H 列 )が 自 動 計 算 される. 謝 意 FFT( 高 速 逆 フーリエ 変 換 )のプログラムは, 露 本 伊 佐 男 氏 のホームページに 記 載 され たものを 使 用 しました.URL は 次 のとおりです. http://tsuyu.cocolog-nifty.com/blog/2007/03/publi.html 逆 フーリエ 変 換 のプログラムは, 上 記 FFT プログラムを 逆 FFT 用 に 修 正 したものです. この 場 を 借 りてお 礼 申 し 上 げます. 参 考 : ボタン 1 には, 以 下 のプログラムを 記 述 する.
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-11/14 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 '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 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
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-12/14 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 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 'データの 出 力 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)) ' 逆 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) '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 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 埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-13/14
埼 玉 工 業 大 学 技 術 資 料 ( 小 西 克 享 ) Excel によるデジタルフィルタ-14/14 j = j - k k = k / 2 Loop j = j + k Cells(i0, 8) = xr(0) / n 'データの 出 力 For i = 2 To n Cells(i0 + n - i + 1, 8) = xr(i - 1) / n End Sub http://www.sit.ac.jp/user/konishi/jpn/tech_inform/pdf/digitalfilterbyexcel.pdf Copyright c 2013-2014 小 西 克 享, All Rights Reserved. 個 人 的 な 学 習 の 目 的 以 外 での 使 用, 転 載, 配 布 等 はできません.