macroscope

( はてなダイアリーから移動しました)

時系列を振動数によってわける -- ディジタルフィルター (とりあえずの覚え書き)

【まだ書きかえます。どこをいつ書きかえたか、かならずしも明示しません。】

- 1 -
数量の時系列を、正弦波 (sin, cos) のかさねあわせ (線形結合) としてとらえることは、おおくの専門でつかわれる方法だ。気象の分野では、時系列をたくさんの振動数の成分からなるものとかんがえて、(単一の振動数の成分をとりだすのではなく) 振動数の区間ごとにわけることがおおい。【[注] 「振動数」と「周波数」は同じことなので「振動数」で統一しておく。】

そのような方法について、わたしは、1980年代前半、気象学の大学院生だったころ、助手(いまならば助教)だった新田 勍 [つよし] さんから手ほどきをうけ、いくつかの文献を読んだり、プログラムを動かしてみたりして勉強した。1980年代後半には教えるがわになった。

そして、2000年に地理学の演習の授業をしたときに、このような教材ページをつくった。

そこに「パワースペクトル」と「ディジタルフィルター」を追加しようと思った痕跡がのこっているが、実際には書かなかった。

スペクトル解析については、その後、このブログに[2011-10-22 気象・気候の問題に使われるスペクトル解析についての序説]を書いた。また、新田さんと林 良一[よしかず]さんによる Fortranプログラムをわたしがいくらか見やすく改造したものが、1987年に東京大学 大型計算機センターのプログラムライブラリに登録されたのだが、そのプログラムをわたしのウェブサイトにも置くことにした。上記のブログ記事からプログラムのページへのリンクを、きょう(2020-01-28)追加した。

- 2 -
ディジタルフィルターについては、わたしの知識は、気象学のうちで経験したたぐいの応用を教えることはできても、ひろく一般に読まれる教材を書くには不じゅうぶんだと思うのだが、ともかく知っていることを書きとめておく。

時系列データを、振動数の区間ごとにわけるには、データをフーリエ変換して (パワーをもとめるのではなく、各振動数の sin 成分の係数と cos 成分の係数、あるいは複素指数関数の係数の実数部と虚数部、を両方とっておいて)、振動数の座標軸のうえで区間を選択し、選択した区間の値を逆フーリエ変換して時系列にもどす、という方法もある。これは、理屈はわかりやすい。ただし、時系列が長くなる(時間ステップ数がおおくなる)と、計算量が重くなる。また、現実の時系列には端があるが、フーリエ変換する際には周期性を仮定するので、端のところのデータをゆがめてしまう。その影響は理屈のうえでは時系列全体におよぶ (実質的には端の付近だけに見られるのがふつうだが)。

振動数の座標軸にもっていかないで、時間軸上の操作だけで、振動数の区間ごとに情報をわけたい。つぎのような需要がある。

  • 振動数が高い(周期が短い)ほうをとりだす ... high-pass
  • 振動数が低い(周期が長い)ほうをとりだす ... low-pass
  • 振動数がふたつの値ではさまれた区間をとりだす ... band-pass

これをやる方法がディジタルフィルターだ。ここでいうディジタルは、あたえられるデータが、時間軸上で、連続時間の関数ではなく、離散的な時点列での値としてあたえられているという意味だと思う。(たぶん、ディジタル計算機で処理されるとか、有限桁数の固定小数点数または浮動小数点数で計算されるとかいう意味ではないだろうと思う。) ディジタルフィルターが対象とするのは、ほとんどのばあい、等間隔の時点列であたえられた数値の列だ。

ディジタルフィルターには、大きくわけて、非再帰型 (non-recursive) のものと再帰型 (recursive)のものがある。【[注] recursive を「巡回型」としている文献もある。】

非再帰型のフィルターでは、時刻 t での出力値は、時刻 t 前後の有限個の入力値の線形結合で表現される。入力を x、出力を y とすれば

y(t) = Σk=k1 .. k2 a(k) x(t+k) ... (1)

のような形になる。 (数式の総和のΣによる表現をかきたかったのだが、このブログ プラットフォームでどう指定すればできるのかわからないので、標準からはずれてもうしわけないが、下限・上限を両方とも下つき添え字のように書いた。) 係数の組 a(k) は、時間差 k だけにより、時間 t の経過によっては変化しない。【この時間差 k は、正負の符号を逆にとってラグ (lag、遅れ) というほうがふつうなのだけれど、ここではこうしておく。】
たとえば、

y(t) = (1/4) x(t-1) + (1/2) x(t) + (1/4) x(t+1) ... (2)

は low-pass フィルタということができ、(1)の一般形と対応させれば k1 = -1 , k2 = +1 である。

非再帰型フィルターでは、入力(の特定の部分)の影響のおよぶ範囲が限定される。(1)の式で、x(t1)の影響がおよぶのは、y(t1-k1) から y(t1-k2)までなのだ。この観点で、非再帰型フィルターを、finite impulse response (FIR、有限インパルス応答) フィルターともいう。

再帰型のフィルターでは、時刻 t での出力値は、時刻 t 前後の有限個の入力値と、時刻 t の前の有限個の出力値の線形結合で表現される。

y(t) = Σk=k1 .. k2 a(k) x(t+k) + Σk=k3 .. -1 b(k) y(t+k) ... (3)

のような形になる。再帰型フィルターでは、入力(の特定の部分)の影響は、(だんだんうすまるが)それ以後の全体に出る。この観点で、再帰型フィルターを、infinite impulse response (IIR、無限インパルス応答) フィルターともいう。

再帰型フィルターは、端の影響が全体におよぶこと、係数の組がどうしても非対称になるので出力の波の位相が入力のその振動数の成分からずれること、などの欠点がある。しかし、係数の項数がすくなくてすみ、端の付近でデータ点数をあまりうしなわないですむという長所がある。なお、位相のずれは、リアルタイムの利用でないばあいは、同じ係数のフィルターを時間前向きとうしろ向きで2回かけることによってなくすことができる(つぎにのべる Murakami 1979はそうしている)。

フィルターの理想は、振動数軸上で、とおしたい区間の成分は全部とおし、とめたい区間の成分は全然とおさないものなのだが、現実のディジタルフィルターでは、いくらかまざってしまう。項数をふやしたほうが、特性は理想にちかづく。しかし端の付近でうしなうデータ点数がふえることや、計算量がふえることなどの欠点もある。したがって、目的によって、いろいろな項数のフィルターがつかわれることになる。

- 3 -
(これから学ぶ人にすすめられるかどうかわからないが) わたしが学んだ文献をあげておく。

1980年ごろ、気象学の分野の研究論文で、ディジタルフィルターを使い、その手つづきの説明をつけている文献として、Blackmon (1976)とMurakami (1979) があった。Blackmon は非再帰型、Murakami (村上 勝人[まさと] )は再帰型のフィルターを使っていた。

Murakami は、地震波による地中構造の探査の分野の論文雑誌にのった Shanks (1967) の論文を参考にしていた。Shanksがのべた一連の方法のうちで、Murakami は、いちばん簡単な、2ステップ前までの値だけを使う方法をとった。なお、Shanks (1967) の話題については、斉藤・石井 (1969) による日本語の解説もある。

わたしは、Murakami (1984) の雲量の季節内変動に関心があったので、そこでもつかわれていたMurakami (1979) のディジタルフィルターを実行する計算プログラムを組んだが、自分ではそれをつかった研究成果を出さなかった。しかし後輩にあたる高薮 縁[ゆかり]さんが修士論文としてとりくんでいた熱帯の大気の波動の研究で同じ方法をつかうのを見てはいた。その修士論文は Nitta & Takayabu (1985) の論文になった。

ディジタルフィルターの基本概念については、わたしは Hamming の本の日本語版で勉強した。また、ディジタルフィルターやパワースペクトルをふくむ、時系列データの統計的処理については、Bendat & Piersol の本の日本語版や、得丸ほか (1982) の本で勉強した。

文献

  • Julius S. Bendat & Allan G. Piersol, 1971: Random Data: Analysis and Measurement Procedures. Wiley.
  • [同、日本語版] J. S. ベンダット、A. G. ピアソル 著, 得丸 英勝 ほか 訳, 1976: ランダムデータの統計的処理。培風館。とくに 9.2節 数値フィルタ法。
  • Maurice L. Blackmon, 1976: A climatological spectral study of the 500 mb geopotential height of the Northern Hemisphere. Journal of the Atmospheric Sciences, 33: 1607-1623. doi: 10.1175/1520-0469(1976)033<1607:ACSSOT>2.0.CO;2
  • R. W. Hamming, 1977: Digital Filters. Prentice-Hall.
  • [同、日本語版] R. W. Hamming 著、宮川 洋、今井 秀樹 訳, 1980: ディジタル・フィルタ。科学技術出版社。
  • Masato Murakami, 1979: Large-scale aspects of deep convective activity over the GATE area. Monthly Weather Review, 107: 994-1013. doi: 10.1175/1520-0493(1979)107<0994:LSAODC>2.0.CO;2 (公開)
  • Masato Murakami, 1984: Analysis of the deep convective activity over the Western Pacific and Southeast Asia. Part II: Seasonal and intraseasonal variations during northern summer. Journal of the Meteorological Society of Japan, 62: 88-108. https://doi.org/10.2151/jmsj1965.62.1_88 (公開)
  • Tsuyoshi Nitta & Yukari Takayabu, 1985: Global analysis of the lower tropospheric disturbances in the tropics during the northern summer of the FGGE Year. Part II: Regional characteristics of the disturbances. Pure & Applied Geophysics (PAGEOPH), 123: 272-292. https://doi.org/10.1007/BF00877023 (公開)
  • 斉藤 正徳, 石井 吉徳, 1969: 簡単なRecursiveフィルター。物理探鉱 22(6): 527-532. http://www.segj.org/~kashima_admin/B_list.html#V22N6 (会員制)
  • J. L. Shanks, 1967: Recursion filters for digital processing. Geophysics, 32: 33-51. https://doi.org/10.1190/1.1439855 (有料) (Murakami 1979 の参考文献リストでは雑誌名が Geophysica となっているが、まちがい。)
  • 得丸 英勝 ほか 著, 1982: 計数・測定 -- ランダムデータ処理の理論と応用 (工学基礎講座 18)。培風館。