フーリエ変換とベイズモデルおよびトレンドモデルによる時系列データのノイズ除去

小野寺 光永, 井須 芳美, 長嶋 雲兵, 吉田 裕亮, 細矢 治夫, 永川 祐三

Return

1 はじめに

　時間の経過とともに不規則に変動する現象の時系列データは、自然現象をはじめ、医学、経済等様々な分野で見ることができる。その解析により我々は多くの情報を得ることができるが、一般に時系列データの解析では、観測された時系列データからノイズを除去して真のデータを引き出すことが重要な問題となる。特にSN比が悪くなるような大きなノイズがのる場合、それが解析の著しい障害となる。例えば生体の断層イメージを測定する核磁気共鳴(NMR)のシグナルや、Belousov-Zhabotinsky反応(BZ)などの振動反応の酸化還元ポテンシャルには、測定環境によりSN比が急激に悪くなる時系列データが見られ、そのノイズ除去が問題となっている。また生体における時系列データの一つである心電図においても、携帯心電計による測定ではSN比が急激に悪くなるようなノイズがのるケースがあり、それが狭心症などを診断する心電図解析の一つであるRR間隔測定の際に大きな障害となる。
　ノイズ除去の方法として、これまでにフーリエ変換[1]、ベイズモデル[2, 3]、トレンドモデル[7]等が提案されてきた。フーリエ変換による方法は、時系列データを正弦波の和として表すことにより周波数分解し、その高周波成分や低周波成分の一部を除くことによりノイズの除去を行う。ベイズモデルでは、観測値に対する確率分布を導くことにより真のモデルの推定を行う。トレンドモデルは、パラメータの確率的変化を表現したモデルであり、時系列の長期的変動であるトレンドを推定する。これらの手法のうちベイズモデルやトレンドモデルは、精度の良い結果を得るために統計的知識を伴った多くの試行錯誤を必要とする。一方フーリエ変換による方法は、統計的な予備知識を必要としない反面、特定の周波数成分を除くことによるトレンドへの影響を見る必要がある。
　そこで本研究では、フーリエ変換、ベイズモデル、トレンドモデルの3手法を用いて、SN比が急激に悪くなるデータを含めた2種類のサンプルデータのノイズ除去を行い、それぞれのモデルによる結果の比較、検討を行った。

2 サンプルデータ

　サンプルデータとして以下の2種類のデータを用いた。それぞれを簡単に説明する。

2. 1 Multiple nonstationary frequencies

　NMRやBZ反応のシグナルは、Multiple nonstationary frequenciesで表されることが知られている。そこで本研究では、Multiple nonstationary frequenciesの近似データを式 1 により作成し、この時系列データにSN比がそれぞれ20, 10, 5となるようノイズを加えた3種類の時系列データをサンプルデータとして用いた。その中のSN比が5のデータをFigure 1に示す。ここで、SN比 sn はデータf(t)の標準偏差 をノイズの標準偏差 pn で割った値で、式 2 のように表すことができる。


 Figure 1. Multiple nonstationary frequencies(SN比=5)

　式 1 により作成されたデータは、Figure 1からわかる通り、時間における振幅の変化が大きく、データ中の時間変化によるSN比の差が大きい。そこでFigure 1の縦線で示したように、データを[0,255], [256,511], [512,767], [768,1023]の4つの時間区間に分け、データ中のSN比の大きさの変化による結果の違いを見ることにした。SN比が20, 10, 5であるそれぞれのノイズデータの時間区間ごとのSN比をTable 1に示す。[0,255]は区間(0)、[256,511]は区間(1)、[512,767]は区間(2)、[768,1023]は区間(3)とした。SN比が5のノイズデータは、区間[3]でSN比がかなり悪くなっており、シグナルよりノイズの方が大きくなっている。

Table 1 各区間のSN比
全体のSN比 区間のSN比
(0)(1)(2)(3)
2060.1530.376.572.69
1030.0715.183.291.35
515.047.591.640.67

2. 2 心電図

実際の時系列データへの応用として、携帯心電計により測定された心電図のデータを使うことにした(Figure 2)。一般の心電図検査では仰臥位安静時に測定を行うが、携帯心電計では、日常生活中の運動負荷の状態変化をみることができ、こうして心電図の監視を継続的に行うことで、高血圧障害などの患者の在宅医療が進むものと期待されている。しかし、携帯心電計により測定されたデータには、携帯心電計から記録装置への電波の伝送ミスのため、Figure 2中のt=100, 360付近で見られるようにSN比が急激に悪くなる大きなノイズが含まれる。このようなノイズは、RR間隔（心電図内の各QRS波の間隔(Figure 3）測定の障害となる。この間隔は運動により変化するが、心臓病との関係も大きく、期外収縮ではその短縮、房室ブロックではその延長が見られる。そのためRR間隔時系列は、医療における重要な診断材料となっている[4]。従って、その正確な測定はたいへん重要であり[5]、そのためのノイズ除去が必要となる。
　またこのような測定時のノイズは、BZ反応などにおいても見ることができるので、心電図以外のデータのノイズ除去への応用分野も広い。


 Figure 2. 携帯心電計により測定された心電図


Figure 3. 心電図のRR間隔

3 ノイズ除去の方法

3. 1 フーリエ変換によるノイズ除去

　フーリエ変換では、時間関数(g(t), t=0,…,N-1)から周波数関数(G(f))を求める（式5）。逆フーリエ変換では、反対に周波数関数から時間関数を求める（式6）。本研究では、フーリエ変換を行い周波数関数を求め、その高周波成分を除くことによりノイズを除去した後、その結果に逆フーリエ変換を行い、時間関数に戻すことを試みた。

　本方法は、ロー／ハイパスフィルタリングとして知られている。

3. 2 ベイズモデルによるノイズ除去

　ベイズモデルによるノイズ除去には田辺らの方法[6]を用いた。以下にその方法を説明する。
　推定すべき関数 f(x)を等間隔な離散点 xj上の値として fj = f(xj)と表現し、fを式 7 のように表す(tは転置)。以下では関数とベクトル fを同一視する。また今、N個のデータ yi(1 <= i <= N)が得られているとし、yiを式 8 のように表す。n は N より十分大きくとる。離散点xj上に常にyiがある必要はない。

　データ yiの誤差εiをεi 〜 i.i.d.N(0,σ2)とすると、fj上のデータである yiの分布は式 9 で表すことができ、この時 y の分布は式 10 で表すことができる。ここで i.i.d.は、independently and identically distributed（定常独立）の意味であり、||・||は、ユークリッドノルムを表す。また E は y と f の対応を表す(N×n)行列で、fit-point xj上にデータ yiがあれば、そのi 行 j 列は1、その他は0となるような行列である。

　次に、f の滑らかさを二階差分の大きさで表現する。各点での二階差分の大きさを式 11の ような量で測り、f の滑らかさの程度を|| dDf ||2で測る。D は式 12 のようにとる。

　ここで、超パラメター d を持つimproperな事前分布式 13 を導入し、f の滑らかさをこの事前分布に含まれる超パラメータ d で制御する。l = n -2、ψは DtD の非零固有値である。

　従って、データ yiが与えられた時の周辺尤度 L(σ2, d )は式 14 のように定義される。

　赤池はベイズ型情報規準量ABIC=-2log L(σ2, d )を最小にするσ2, d を f の推定に用いる方法を提案している。今、ABICを最小にするσ2, dを （最尤推定値）とすると、ベイズの定理により f の事後分布は式 15 で表されるので、この式を最大にする f を最適近似関数 として選ぶことにする。

　ここで p(f|y)は、y という事象が起こった後にfが起こる確率を意味する。L (σ2, d )は、α = d σと超パラメターの変数変換を行うと式 16 のように表すことができる。ここで、|| b- Zα f* ||2は最小二乗問題 minf || b - Zα f ||2の最小二乗誤差で、f*はその時の最小二乗解である。つまり、 は のときのf*に他ならず、問題は式 16 を最大にするσ2, αを求めることになる。

　より式 17 を得ることから、式 16 は式 18 となる。従って、αについてその尤度方程式を解けば良いが、これは解析的に解くことが困難である。そこでαを数値的に求める。いろいろなαに対して式 19 を計算し、そのABICが最小になるようなαを求め、そのαに対応するf*を推定関数にする。

3. 3 トレンドモデルによるノイズ除去

　トレンドは時系列のおおよその傾向を表す。逆にいえば、実際に観測される時系列はトレンドに様々な変動成分が加わったものである。その中で最も単純な場合は、wnを平均0、分散σ2の正規白色雑音として、式 20 のように表現できる。

　トレンド tnは説明変数 xnの多項式で、式 21 のように表されるものとする。

　ここで、式 21 の多項式をより柔軟な関数に拡張する。一般に、k-1次の多項式は k 階の差分方程式(式22)の解とみなすことができるが、差分方程式の代わりに Δktn0 が成り立つように k 階の確率差分方程式(式23)を解くことにする。ただし、Δ≡1-B は、Δtn = tn - tn-1によって定義される時間差分オペレータとし、vn は、平均 0、分散τ2の正規白色雑音とする。式 23 をトレンド成分モデルと呼ぶ。

　従って、時系列 ynからトレンド tnを推定するためには、式 23 と式 20 を組み合わせた次のようなモデルを考えれば良い。これをトレンドモデル(trend model)と呼ぶ。

　式 25 は、時系列 ynをトレンドに独立なノイズが加わったものを観測した結果とみなした観測モデルである。一方式 24 は、そのトレンドの変化の仕方をモデル化したトレンド成分モデルである。
　一般に、トレンド成分モデル式 24 は、式 26 となるので、ci = (-1)i+1kCiによって係数 ciを定義すると、式 27 のように表現できる。

　このモデルは定常性を満たさないが形式的には k 次のARモデルとみなせるので、

　とおくことによりトレンド成分モデルの状態空間表現

　が得られる。これを利用することにより、トレンドモデルの状態空間表現(式30)が得られる。

　従って、トレンドモデルの次数 k 及び分散τ2, σ2が定まると、カルマンフィルタ及び平滑化のアルゴリズムにより状態ベクトルの平滑値 x1|N , … ,xN|Nを求めることができる[7]。また、σ2はカルマンフィルタを用いることによりその最尤推定値が得られることから、本研究では、次数kと分散τ2のみをパラメータとして選択した。

4 結果

4. 1 Multiple nonstationary frequencies

　2.1節で説明したMultiple nonstationary frequenciesの近似関数のノイズ除去をフーリエ変換、ベイズモデル、トレンドモデルそれぞれにより行った。各手法のパラメータの数値は、それぞれのSN比のノイズデータで区間(0)の相対誤差が最小になるように選択した。つまりフーリエ変換では、サンプル数として2048点を用い、式 5 , 6の k = 0〜410の G(k)を0と置き、高周波数成分を除いた。ベイズモデルでは、3.2節で説明した滑らかさを表すパラメータαをデータのSN比が20のときα=7、10のときα=12、5のときα=16とした。トレンドモデルでは、3.3節で説明したモデルの次数を1、トレンドのノイズの分散τ2をデータのSN比が20のときτ2=1.9、10のときτ2=1.6、5のときτ2=0.9とした。
　結果をTable 2に示す。Table 2は左列から順に、全体のSN比、(0)〜(3)の各区間とそのSN比、式 1 に示したノイズを加える前のデータとノイズを加えた後のデータとの各区間における相対誤差の平均値(|e-s|/e)、ノイズを加える前のデータとフーリエ変換を用いてノイズを除去したデータとの各区間における相対誤差の平均値(|e-f|/e)、同様にベイズモデルによる結果の相対誤差の平均値(|e-b|/e)とトレンドモデルによる結果の相対誤差の平均値(|e-t|/e)をそれぞれ表す。また、Figures 4 - 6にSN比が5のデータをフーリエ変換、ベイズモデル、トレンドモデルそれぞれによりノイズ除去した図を示す。Figure 7は、区間(2)における拡大図である。

Table 2 フーリエ変換(f)、ベイズモデル(b)、トレンドモデル(t)の各区間における実データとの相対誤差
全体のSN比区間区間のSN比|e-s|/e|e-f|/e|e-b|/e|e-t|/e
20(0)60.150.100.040.040.08
(1)30.370.120.040.050.08
(2)6.574.173.173.053.34
(3)2.691.511.731.001.18
10(0)30.070.190.070.080.13
(1)15.180.250.080.080.13
(2)3.298.346.175.556.55
(3)1.353.012.571.842.30
5(0)15.040.380.130.140.22
(1)7.590.490.150.150.22
(2)1.6416.6912.1710.4012.01
(3)0.676.034.293.544.22

Table 3 フーリエ変換(f)、ベイズモデル(b)、トレンドモデル(t)の全区間における実データとの相対誤差
全体のSN比|e-f|/e|e-b|/e|e-t|/e
201.241.041.17
102.221.892.28
54.193.564.17

Table 2の右から4列目に示したノイズを加える前のデータと加えた後のデータとの相対誤差(|e-s|/e)と、各手法による相対誤差を比較することで、各手法によるノイズ除去の程度を見ることができる。|e-s|/eと比較して相対誤差が小さいほど、ノイズが除去されたことを表し、逆に大きくなっていればノイズの除去に失敗したことを表す。
　まずフーリエ変換の相対誤差は、区間(0), (1)ではSN比=20, 10, 5全てのデータで|e-s|/eの1/2以下になっている。区間(2), (3)では区間(0), (1)ほどではないが、SN比=20の区間(3)を除いて相対誤差が減少している。
　ベイズモデルの相対誤差は、区間(0), (1)では、フーリエ変換の結果と同様に全てのデータで|e-s|/eの1/2以下になっている。また区間(2), (3)では、区間(0), (1)ほどではないが、全てのデータで相対誤差が減少している。
　次にトレンドモデルの相対誤差を見ると、全ての区間で相対誤差が減少しているが、区間(0), (1)の相対誤差がフーリエ変換やベイズモデルと比較すると大きい。
　各手法の相対誤差を比較すると、区間(0), (1)ではフーリエ変換とベイズモデルがほぼ同じ誤差であり、トレンドモデルと比較して誤差が小さい。区間(2), (3)ではベイズモデルの誤差が最も小さい。また、Table 3の全区間における誤差を比較してもベイズモデルの相対誤差が最も小さい。
　Figures 4 - 6の時刻512以降(区間(2),(3))をみると、Figure 6のトレンドモデルによる結果はFigures 4, 5に示した他の手法による結果と比較して滑らかさに欠けている。これは区間(2)の拡大図であるFigure 7をみると一層よくわかる。τ2を小さい値にすると滑らかな結果が得られるが、本研究では、τ2を区間(0)の相対誤差により選んだのでこのような結果が得られた。一方、同様にパラメータ(σ2)を選択したベイズモデルは、Figure 7に示すように、3手法の中で最も滑らかな結果が得られていることがわかる。
　これらの結果から、このサンプルデータによる実験では、ベイズモデルが最もSN比によらず精度の良いノイズ除去を行うことができるという結論を得た。


 Figure 4. ノイズデータとフーリエ変換(fft)によるノイズ除去(SN比=5)


Figure 5. ノイズデータとベイズモデル(bayes)によるノイズ除去(SN比=5)


Figure 6. ノイズデータとトレンドモデル(trend)によるノイズ除去(SN比=5)


Figure 7. フーリエ変換(fft)、ベイズモデル(bayes)、トレンドモデル(trend)によるノイズ除去(SN比=5)

4. 2 心電図

　Figure 2に示した心電図を用いて数値実験を行った。この心電図は、健康な人の階段昇降時の心拍状態を測定したものであり、データ数は512である。
　まずt=100, 360付近の電波の伝送ミスによるノイズを除くため、フーリエ変換によるノイズ除去を行った。その結果をFigure 8に示す。ここでは、サンプル数512で、式 5, 6の k =0〜170の G(k)を0と置き、高周波数成分を除いた。QRS波のピークがFigure 2と比較してほとんど変化がないのに対し、t=100, 360付近で見られた大きな振れは抑えられていることがわかる。
　また、Figure 2では運動負荷による基線の揺れが見られるが、このこともRR間隔時系列を取り出す際に障害となる。そこで、基線を一定にするため、高周波成分に加え、低周波成分を除くことを試みた。その結果をFigure 9に示す。ここでは、式 5, 6の k =508〜512の G(k)を0と置き、低周波数成分を除いた。Figure 9を見ると、Figures 2, 8と比較して基線の揺れが抑えられ、一定となっていることがわかる。
　次に、Figure 2をベイズモデルにより平滑化した結果をFigures 10, 11に示す。Figure 10はα=1.0、Figure 11はα=5.0である。αの値が大きいほど滑らかな曲線が得られるため、Figure 10ではt=100, 360付近の大きな振れが残っているが、Figure 11ではほぼ完全に除かれている。またFigure 10とFigure 11のQRS波を比較すると、Figure 11の方が多少波が低くなっているが、QRS波を捉えるには十分である。
　次に、Figure 2をトレンドモデルにより平滑化した結果をFigures 12, 13に示す。Figure 12はτ2=1.4、Figure 13はτ2=0.5である。τ2の値が小さいほど滑らかな曲線が得られるため、Figure 12ではt=100, 360付近に振れが残っているが、Figure 13ではほぼ完全に除かれている。
　Figures 9, 11, 13を比較すると、Figure 9は全体に細かい波がのっており、そのため滑らかさに欠けるが、このような細かい波は、RR間隔測定の障害にはならない。またFigure 9で示した通り、フーリエ変換による方法では、基線の揺れも抑えられるため、Figure 9にのっている細かい波が気になるならば、ベイズモデルやトレンドモデルにより平滑化をした後、フーリエ変換を用いて基線を一定にする方法が有効であると考えられる。


 Figure 8. フーリエ変換によるノイズ除去 (高周波成分を除去)


Figure 9. フーリエ変換によるノイズ除去(高、低周波成分を除去)


Figure 10. ベイズモデルによるノイズ除去(α=1.0)


Figure 11. ベイズモデルによるノイズ除去(α=5.0)


Figure 12. トレンドモデルによるノイズ除去(τ2=1.4)


Figure 13. トレンドモデルによるノイズ除去(τ2=0.5)

5 まとめと今後の課題

　本研究では、2種類のサンプルデータのノイズ除去をフーリエ変換、ベイズモデル、トレンドモデルの3手法により行い、各手法による結果の違いを検討した。
　まずNMRやBZ反応のシグナルに見られるMultiple nonstationary frequenciesの近似関数にSN比が20, 10, 5となるようノイズを加え、そのノイズ除去を行った。このデータは時間によるSN比の変化が大きいため、データを4つの時間区間に分け、それぞれの区間の誤差を見た。その結果、SN比が大きいはじめの2つの区間では、フーリエ変換とベイズモデルにより、ノイズを除去した後の相対誤差をノイズを除去する以前の相対誤差の1/2以下にすることができた。トレンドモデルは、フーリエ変換及びベイズモデルと比較して誤差が大きかった。SN比が小さい残りの2区間では、ベイズモデルが最も誤差が小さかった。また、全区間の誤差を比較してもベイズモデルによる結果が最も誤差が小さかった。これらの結果から、このサンプルデータによる実験では、ベイズモデルが最もSN比によらず精度の良いノイズ除去を行うことができるという結論を得た。
　次に、携帯心電計により測定された心電図を用い、心電計と記録装置間の電波の伝送ミスによるノイズ除去を行った。その結果、フーリエ変換、ベイズモデル、トレンドモデル、どの手法を用いても、電波の伝送ミスによる、SN比が大きいノイズを除くことができた。さらにフーリエ変換では、運動負荷による基線の揺れも抑えることができた。しかし、フーリエ変換による方法では、全体に細かい波がのり滑らかさに欠けるため、ベイズモデルやトレンドモデルにより平滑化をした後、フーリエ変換を用いて基線を一定にする方法が有効であると考えられる。これらの手法を用いたノイズ除去により、医療における重要な診断材料である心電図のRR間隔の測定が容易になることが期待される。
　ノイズ除去の手法ごとの長所、短所が、Multiple nonstationary frequenciesの近似関数や、本来の心電図などの「関数群」によるものなのか、それともその上にかぶせた、ランダムな雑音や、電波伝送ミスにより起こる大きなノイズ等の「ノイズの種類」によるものなのか、という疑問が残るが、この究明を今後の課題として挙げたい。また、心電図以外の時系列データにおけるノイズ除去を行い、その効果を確かめていきたい。

参考文献

[ 1] 安居院猛, 中嶋正之, FFTの使い方, 秋葉書店, 東京 (1994).
[ 2] 鈴木雪夫, 国友直人編, ベイズ統計学とその応用, 東京大学出版会, 東京 (1992).
[ 3] 繁桝算男, ベイズ統計入門, 東京大学出版会, 東京 (1994).
[ 4] 五島雄一郎, 大林完二監修, 心電図のABC, 日本医師会, 東京 (1997).
[ 5] 早野順一郎, 心拍変動による自律神経機能解析, 循環器疾患と自律神経機能, 医学書院, 東京 (1996), p58.
[ 6] 田辺國士, 田中輝雄, ベイズモデルによる曲線・曲面のあてはめ, 月刊地球, 5, 179-186 (1983).
[ 7] 北川源四郎, 時系列解析プログラミング, 岩波書店, 東京 (1993).
Return