気象学でハイパスにFFTの周波数カットでなく移動平均の差し引きを使う理由

問いは木下が出したもので、答えは Claude Code と調べて書いたものです。

昨日気象学の研究してD進した人と飲んでて、気象ではよく移動平均を引くことでハイパスにして短周期の特性をみつけるんだよねっていっていたんだけど、なんでFFTとかで周波数カッティングして戻すんじゃなくて移動平均を引くみたいな原始的な方法でやっているの？

答え：FFTで周波数をすっぱり切るほうが、時間領域では「原始的な」フィルタより扱いにくいから。観測データは欠けていて、周期的でもなく、どこにでも急な出来事が入っている

周波数を鋭く切ると、時間方向には長く広がる。理想的なハイパス（ある周波数より下をゼロにする）は、時間領域では $\sin x / x$ が無限に続く重みとの畳み込みと同じである。そのため、台風の通過や寒冷前線のような急な変化があると、その前後に実在しない振動（リンギング、ギブズ現象）がにじみ出る。短周期の現象を探すのが目的なのに、偽物の短周期の波を作ってしまう。

FFT にはほかにも弱点が三つある。①データが周期的だとみなされる。 FFT は記録の終わりと始まりがつながっていると仮定するので、両端の値が食い違うと、そこに大きな段差があるものとして扱われ、両端が汚れる。②欠測を扱えない。 観測所の故障や衛星の欠けでデータには穴があるが、FFT は等間隔で欠けのない系列でないと使えない。移動平均なら、窓の中にある点だけで平均を取れば済む。③結果が局所的でない。 FFT でカットした値は、記録全体のデータの影響を受ける。何十年も離れた年の値が入れ替わるだけで今日の結果が変わり、新しい観測を足すたびに過去の結果まで変わる。移動平均なら、時刻 $t$ の結果は前後 $T$ の範囲だけで決まる。

もう一つの理由は、説明のしやすさと慣習である。「31日移動平均からのずれ」と書けば誰でも同じ計算を再現でき、論文どうしで比べやすい。

ただし、気象学で使うのは素の移動平均だけではない。素の移動平均の伝達関数は $1-\sin(\omega T)/(\omega T)$ で、サイドローブのせいで、ある周波数の成分を約22%持ち上げてしまう。そこで実際の研究では、移動平均の重みを工夫した Lanczos フィルタ（Duchon, 1979）や Butterworth フィルタが定番になっている。どれも「有限の幅の窓で重み付き平均を取る」形は同じで、FFT の矩形カットは窓の幅を無限にした極端な場合にあたる。窓が狭いほど時間には局所的で、周波数の切れ味は鈍くなる。その中間の妥協点を選んでいる、というのが正確な見方である。

補足：周波数を鋭く切るとなぜ sinc の畳み込みになり、なぜ偽の振動が出るのか

「周波数を鋭く切ることは、時間領域では sin x / x が無限に続く重みとの畳み込みと同じです。そのため台風の通過や前線のような急な変化があると、その前後に実在しない振動（リンギング）が出ます。」みたいなのってなんで？

sinc になる理由。 周波数領域で掛け算することは、時間領域では畳み込みにあたる（畳み込み定理）。周波数を切ることは、$|\omega|<\omega_c$ で 1、それより外で 0 になる矩形 $H(\omega)$ を掛けることなので、時間領域では $H$ を逆フーリエ変換した $h(t)$ との畳み込みになる。計算すると $h(t)=\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{i\omega t}\,d\omega=\frac{\sin(\omega_c t)}{\pi t}$ となる。これは、時間の矩形窓（移動平均）を変換すると sinc になる計算を、時間と周波数を入れ替えてやり直したものである。$h(t)$ は、周波数 $0$ から $\omega_c$ までの cos 波を同じ強さで足し合わせたものと言い換えられる。$t=0$ では全部の山がそろって大きなピークになり、離れると位相がばらけて打ち消し合う。ところが、いちばん高い周波数 $\omega_c$ の付近には、それを打ち消す相手（$\omega_c$ より上の成分）がもう無いので、周波数 $\omega_c$ の揺れが $1/t$ でしか減らずに残り続ける。矩形の縁が切り立っているほど、もう一方の領域ではゆっくりとしか減らない。

偽の振動が出る理由。 $h(t)$ には負の重みの谷が交互に並ぶので、これとの畳み込みは「負の重みも混ざった加重平均」になる。一定の値から急に跳ね上がる段差に掛けると、窓の負の谷が段差をまたぐ位置で、平均が本来の値を通り越してしまう。結果は $\frac12+\frac{1}{\pi}\mathrm{Si}(\omega_c t)$（$\mathrm{Si}$ は正弦積分）で、段差の高さの約 9% だけ行き過ぎ、その前後で周波数 $\omega_c$ の振動がゆっくり減衰していく。これがギブズ現象で、矩形波のフーリエ級数を途中の項で打ち切ったときの角の飛び出しと同じものである。周波数で見れば、段差はあらゆる周波数を含み、高い周波数の成分は振動を打ち消して角を鋭くする役を担っている。そこを切り落とすと、打ち消し役がいなくなった振動だけが残る。ハイパスは元のデータからこのローパスを引いたものなので、同じ振動が段差の前後に現れる。

移動平均と比べると。 移動平均の重みはすべて正（矩形）なので、段差に掛けても値はなだらかな坂になるだけで、行き過ぎも振動も起きない。そのかわり、周波数の側にサイドローブを抱える（伝達関数が約22%持ち上がる）。時間で切り立った窓は周波数で揺れ、周波数で切り立った窓は時間で揺れる。どちらかの揺れは必ず残るので、どちらで揺れを許すかの選択になる。

補足2：持ち上がりは周期 21.7 日（31日移動平均のとき）。直すなら FFT の矩形カットでなく、窓をなめらかにする

1.217 とかになってて、移動平均を引くやりかただとここらへんの周波数帯が21.7%くらい強められちゃうので、なんかもうすこし伝達関数いいかんじのやつにしてハイパスしてみてもいいんじゃないかな〜という所感ですね。今時FFTとinv.FFTで離散データのフーリエ変換操作とかすぐ実装できるでしょうし。

$1-\mathrm{sinc}(\omega T)$ の最大は、$\tan x = x$ の最初の解 $x=\omega T\approx 4.493$ で $1.217$ になる。周期にすると窓の幅 $2T$ の約 0.70 倍で、31日移動平均なら周期 21.7 日の成分が約 22% 増幅される。$\pi<\omega T<2\pi$（周期が $T$ から $2T$ の帯）はすべて 1 を超える。

ただし、直す方向は「FFT で矩形に切る」ではない。矩形に切ると、今度は時間の側でギブズ現象（約 9% の行き過ぎと振動）が出るので、揺れの場所を周波数から時間へ移すだけになる。直すべきは窓の形で、窓の重みを全部正かつなめらかにすれば、両方の揺れを小さくできる。いちばん手軽なのは、移動平均を2回かけることである。矩形と矩形の畳み込みは三角形の窓になり、伝達関数は $\mathrm{sinc}^2(\omega T)$ なので、引いたものは $1-\mathrm{sinc}^2(\omega T)$ となり、$0\le\mathrm{sinc}^2\le1$ だから 1 を超えない（数値でも最大 1.000）。3回かけるとガウス窓に近づき、ガウス窓なら伝達関数 $1-e^{-\omega^2\sigma^2/2}$ は単調で、行き過ぎがどこにも無い。経済統計の「2×12 移動平均」や気象の Lanczos フィルタも、この窓の整形の仲間である。FFT で実装するなら、$H(\omega)$ を矩形でなくなめらかな形（たとえば $1-e^{-\omega^2\sigma^2/2}$）にすればよく、それは時間領域でなめらかな短い窓との畳み込みと同じものになる。道具を FFT にするか畳み込みにするかより、伝達関数の形を選ぶことが本題である。端の扱い（折り返し・欠測）は FFT のほうが手間がかかる。