気象情報解析特論　第4回　ローパスフィルタとハイパスフィルタ

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math

下準備：東京の気温

今回の話は，まず東京の気温から始めます。

東京の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv （出典は環境情報論第3回に記述）

In [2]:
tokyo_temp = np.genfromtxt("Tokyo_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
tokyo = tokyo_temp[:, 2]

# 今回は1990年から2019年の30年分のデータを用いる
tokyo = tokyo[(1990 <= y)*(y <= 2019)]
m = m[(1990 <= y)*(y <= 2019)]
y = y[(1990 <= y)*(y <= 2019)]

# 気候値の計算
tokyoc= np.zeros((12))
for mm in range(1, 13): 
    tokyoc[mm-1] = np.nanmean(tokyo[m==mm], 0)

# 偏差の計算
tokyoa = np.zeros((tokyo.shape))
for yy in range(1990, 2020):
    for mm in range(1, 13):
        tokyoa[(y==yy)*(m==mm)] = tokyo[(y==yy)*(m==mm)] - tokyoc[mm-1]
        
# スペクトル解析するときにはデトレンドしましょう
tokyoa = signal.detrend(tokyoa)

# 月別の時系列を2つ描画する関数
def plot_2_mon_time(time_series1, time_series2, lower = -3, upper = 3, init_year=1990, fin_year=2019):
    mon = np.arange(init_year, fin_year+1, 1/12)
    plt.plot(mon, time_series1)
    plt.plot(mon, time_series2, 'r')
    plt.plot(mon, 0*time_series1, 'k')
    plt.xlim(init_year, fin_year)
    plt.ylim(lower, upper)
    plt.show()

# 移動平均を取る関数   
def running_mean(time_series, wn):
    b = np.ones(wn)/wn
    time_series_r = np.convolve(time_series, b, mode="same")
    n_conv = math.ceil(wn/2)
    time_series_r[0] = time_series_r[0] * wn/n_conv
    for i in range(1, n_conv):
        time_series_r[i] = time_series_r[i] * wn/(i+n_conv)
        time_series_r[-i] = time_series_r[-i] * wn/(i + n_conv - (wn % 2))
    return time_series_r

移動平均とローパスフィルタ

前回，「移動平均」という概念を勉強しました。移動平均では，「ゆったりとした変動」を取り出すのでした。

In [3]:
tokyoa_r5 = running_mean(tokyoa, 5)
plot_2_mon_time(tokyoa, tokyoa_r5)#元の時系列

ここで勘の良い方は，「ゆったりとした変動」を取り出すということは，パワースペクトルは低周波側（長周期側）だけ残るのでは？と考えたかもしれません。実際にそうなっているか，確かめてみましょう。

In [4]:
# 10ヶ月，15ヶ月移動平均も計算する
tokyoa_r10 = running_mean(tokyoa, 10)
tokyoa_r15 = running_mean(tokyoa, 15)

# パワースペクトルを計算する関数
def calc_power(time_series, delta_t = 1, ave_num = 10):
    N = time_series.shape[0]
    frequency = np.fft.fftfreq(N, delta_t)[0:math.floor(N/2)]
    time_series_k = np.fft.fft(time_series)[0:math.floor(N/2)] 
    power = np.abs(time_series_k)**2/N;
    
    frequency_mean = np.zeros(int(math.floor(N/2)/ave_num))
    power_mean = np.zeros(int(math.floor(N/2)/ave_num))

    for nn in range (0, int(math.floor(N/2)/ave_num)):
        frequency_mean[nn]= np.mean(frequency[ave_num*nn:ave_num*(nn+1)-1]);
        power_mean[nn] = np.mean(power[ave_num*nn:ave_num*(nn+1)-1]);
    
    return frequency_mean, power_mean

# パワースペクトルを計算
freq_mean_tokyoa, power_mean_tokyoa = calc_power(tokyoa)
freq_mean_tokyoa_r5, power_mean_tokyoa_r5 = calc_power(tokyoa_r5)
freq_mean_tokyoa_r10, power_mean_tokyoa_r10 = calc_power(tokyoa_r10)
freq_mean_tokyoa_r15, power_mean_tokyoa_r15 = calc_power(tokyoa_r15)

plt.plot(freq_mean_tokyoa, power_mean_tokyoa, 'b', label="Without Running Mean") # 移動平均前の東京の気温偏差
plt.plot(freq_mean_tokyoa_r5, power_mean_tokyoa_r5, 'b--', label="5 month Running Mean") # 5ヶ月移動平均後の東京の気温偏差
plt.plot(freq_mean_tokyoa_r10, power_mean_tokyoa_r10, 'r--', label="10 month Running Mean") # 10ヶ月移動平均後の東京の気温偏差
plt.plot(freq_mean_tokyoa_r15, power_mean_tokyoa_r15, 'r', label="15 month Running Mean") # 15ヶ月移動平均後の東京の気温偏差
plt.xlabel("Frequency (/month)") 
plt.ylabel("Power ($°C^2 \mathrm{month}$)")

plt.legend()
plt.show()

たしかに，移動平均のwindowが長くなればなるほど，高周波側が削られて，低周波のみが残っている様子が見えます！

このように，高周波側を削り，低周波のみを通すような操作のことを ローパスフィルタ (low-pass filter) といいます。

移動平均は，もっとも実装の簡単な典型的なローパスフィルタのうちの一つです。

移動平均の除去とハイパスフィルタ

ローパスフィルタがあるなら，高周波のみを通すような操作も考えたくなりますよね。そのような操作のことを，ハイパスフィルタ (high-pass filter) といいます。つまり，背景のゆったりした変動には興味がなくて，時間スケールの短めな変動のみが見たい場合のフィルタリングです。皆さんおなじみのデトレンドも，一応区分としてはハイパスフィルタと言えそうです。

最も簡単なハイパスフィルタは，元の時系列から移動平均を除去することです。すなわちローパスで通ったものだけ通さないようにすると，ハイパスになります。

たとえば，10ヶ月移動平均の時系列を元の時系列から引いてみましょう。

In [5]:
# ローパスを除去（ハイパス）
tokyoa_r10_removed = tokyoa - tokyoa_r10
# 描画
plot_2_mon_time(tokyoa, tokyoa_r10_removed)

こちらは，バタバタしたシグナルだけが残っているので，ちょっと重ねても分かりづらいですね。パワースペクトルを計算してみましょう。

In [6]:
# パワースペクトルを計算
freq_mean_tokyoa, power_mean_tokyoa = calc_power(tokyoa)
freq_mean_tokyoa_r10_removed, power_mean_tokyoa_r10_removed = calc_power(tokyoa_r10_removed)

plt.plot(freq_mean_tokyoa, power_mean_tokyoa, 'b', label="Without Running Mean") # 移動平均前の東京の気温偏差
plt.plot(freq_mean_tokyoa_r10_removed, power_mean_tokyoa_r10_removed, 'r', label="10 month Running Mean") # 10ヶ月移動平均後の東京の気温偏差
plt.xlabel("Frequency (/month)") 
plt.ylabel("Power ($°C^2 \mathrm{month}$)")

plt.legend()
plt.show()

確かに，低周波のところだけちょうど切られています。ハイパスフィルタリングに成功しているようです。

フィルタの応答関数

ローパスフィルタには，実は移動平均以外にも色々な種類があります。というより，そもそも移動平均だけでも，色々なwindowの長さがあります。

そこで，構成したフィルタの特性を示す方法として，応答関数(response function) を描画するという方法があります。応答関数とは，「元の時系列のどの周波数帯でどのくらいフィルタがシグナルを増幅または減衰させるか」を表す関数です。

ある周波数において，応答が1だとそのまま通している，0.5だと半分通している，1.2だと20%増幅している，といった具合です。

たとえば，5ヶ月，10ヶ月，15ヶ月移動平均フィルタの応答関数は次のようになります。

In [7]:
delta_t = 1 # 月別データのサンプリング間隔は1ヶ月
for wn in range(5, 20, 5): # 移動平均の数
    taps = np.full(wn, fill_value=1.0/wn) # 平均するときの重み係数（単純な移動平均では重みは全て等しい）
    w, h = signal.freqz(taps) # 応答関数の横軸の値（w=オメガ; 角周波数）と縦軸の値（h; 応答の大きさ）を計算する関数
    plt.plot(w*1/(2*np.pi*delta_t), np.real(h), label="%d month running mean" % wn)
    plt.legend()
plt.xlabel("Frequency (/month)") 
plt.ylabel("Response Function")
plt.show()

移動平均のwindowが広くなるほど，高周波側が厳しく減衰させられている（より周波数の小さなシグナルのみが通っている）ということがわかります。言い換えると，移動平均のwindowが大きくなるほど，強いローパスフィルタであることを示しています。

一方，気になる点もあります。それは，ローパスフィルタのはずなのに，高周波の方が完全にはゼロになっておらず，バタバタしているということです。これは，移動平均は完全無欠なローパスフィルタではないということの表れです。

一般に，離散的で有限な時系列に対して，完全無欠なローパスフィルタやハイパスフィルタを作るのは困難です。そこで，先人が様々な工夫を凝らして，色々なフィルタを考案してきました。気象データ解析で神山がよく見るのは，ランチョスフィルタ（Lanczos filter) とバタワースフィルタ(Butterworth filter) の二つです。今回は，バタワースフィルタのみ試してみます。

バタワースフィルタ

Pythonでは，scipyのsignalモジュールから，signal.butterでとても簡単にバタワースフィルタを作ることができます。

ここで出てくるカットオフ周波数とは，どの周波数よりも低い（高い）周波数を通すローパス（ハイパス）フィルタを作るかという値です。

また，フィルタオーダーとは，「ある時刻の値を決めるときに周辺時刻のデータをいくつ使うか」です。オーダーが大きいほど，シャープなカットオフが実現されますが，端点では情報が足りないので挙動がおかしくなります。フィルタオーダーが$n$のとき，「$n$次のバタワースフィルタ」のように呼びます。

発展的な注：以下に出てくるbとかaが何を表しているか，そもそもバタワースフィルタはどのような仕組みで動いているのか，などを知りたい方は，「無限インパルス応答フィルタ（IIR filter）」について勉強してください。ただし，学部3年程度の複素関数論の知識が必要になりますので，未習の方は結構長い道のりです。気象データ解析としては，「なぜバタワースフィルタを使うと嬉しいか」「どうやって気象データに当てはめるか」を理解してくれれば十分だと考えています。

In [8]:
# 引数は，時系列，サンプリング間隔（月別データだと1ヶ月），カットオフ周波数，フィルタオーダー
def butterworth_lowpass(time_series, delta_t, f_cut, order):
    b, a = signal.butter(order, f_cut, btype='low', fs = 1/delta_t) # フィルタオーダーは1
    time_series_lp = signal.filtfilt(b, a, time_series) # 与えられた時系列に大してフィルターを実行する
    # 発展的な注：filtfiltは順方向と逆方向の両方でフィルタすることで位相のずれないフィルタにする関数
    return time_series_lp

では，早速東京の気温偏差に，1次のローパスバタワースフィルタをかけてみましょう（カットオフ周波数0.08/month）。

In [9]:
tokyoa_bw_low = butterworth_lowpass(tokyoa, 1, 0.08, 1)
plot_2_mon_time(tokyoa, tokyoa_bw_low)

バタワースフィルタをかけた方が，バタバタしていません（？）。うまくローパスできていそうです。

それでは，肝心の応答関数を見てみましょう。

In [10]:
delta_t = 1

# 比較のため，まず青で10ヶ月移動平均の応答関数を描画
taps = np.full(5, fill_value=1.0/5) # 平均するときの重み係数（単純な移動平均では重みは全て等しい）
w, h = signal.freqz(taps)
plt.plot(w*1/(2*np.pi*delta_t), np.real(h), 'b:', label="5 month Running Mean")

# 赤で1次バタワースフィルタの応答関数を描画
b, a = signal.butter(1, 0.08, btype='low', fs = 1/delta_t)
w, h = signal.freqz(b, a)
plt.plot(w*1/(2*np.pi*delta_t), np.real(h), 'r', label="1st-Order Low-pass Butterworth Filter") 

plt.legend()
plt.xlabel("Frequency (/month)") 
plt.ylabel("Response Function")
plt.show()

移動平均よりもはるかに美しいローパスのカットオフが実現されています！

課題4（締切：次回の授業開始前まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. ローパスフィルタとハイパスフィルタについて，その役割を説明してください。

B問題. バタワースフィルタのカットオフ周波数を色々と変えてみると，ローパスした東京気温偏差の時系列とパワースペクトルはどのように変わるか確かめてください。また，フィルタの応答関数はどのように変わるか確かめてください。

C問題. 1次のハイパスバタワースフィルタを構成し，ハイパスした東京気温偏差の時系列とパワースペクトル，フィルタの応答関数を描いてください。

ヒント：btype='low'をbtype='high'に変えて，あとは資料通りにやるだけです。

D問題. バタワースフィルタ・ベッセルフィルタ・チェビシェフフィルタ・エリプティックフィルタの4つを構成し，それぞれの応答関数を描いてください。また，応答特性の違いを考察してください。

ヒント：https://docs.scipy.org/doc/scipy/reference/signal.html

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論4，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。