import numpy as np
import matplotlib.pyplot as plt
from scipy.special import factorial
import math
from scipy import signal
下準備：SSTデータセットの読み込みと関数の定義¶
今回，Niño3.4指数を使用します。なんのことだか忘れてしまった方はhttp://web.is.ocha.ac.jp/~tsubasa/env_info/Lecture6.html で復習してから臨んでください。
まず，使用データ https://web.is.ocha.ac.jp/~tsubasa/env_info/detrended_ssta_OISST.npz をダウンロードして，この.ipynbファイルと同じディレクトリに置いてください。
loadfile = 'detrended_ssta_OISST.npz' # デトレンド済みの海面水温偏差の入力ファイル名を定義
ssta_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します
ssta = ssta_dataset['ssta'] # 海面水温偏差を変数sstaに保存
lon2 = ssta_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2 = ssta_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = ssta_dataset['y'] # 年(year)を変数yに保存
m = ssta_dataset['m'] # 月(month)を変数mに保存
# 1月始まり12月終わりになるように，1982年から2019年のみのデータを用いる
ssta = ssta[:, :, (1982 <= y)*(y <= 2019)]
m = m[(1982 <= y)*(y <= 2019)]
y = y[(1982 <= y)*(y <= 2019)] # y自身のサイズが変わってしまうので，yは一番最後に書き換えないとダメ
# 領域平均をとる関数
def aave(west, east, south, north, var = ssta):
var = var[(west<=lon2[:, 1])*(lon2[:, 1]<=east),:, :]
var = var[:, (south<=lat2[1, :])*(lat2[1, :]<=north), :]
aave_var = np.nanmean(np.nanmean(var, 0), 0)
return aave_var
# 月別の時系列を描画する関数
def plot_mon_time(time_series, lower = -3, upper = 3, init_year=1982, fin_year=2019):
mon = np.arange(init_year, fin_year+1, 1/12)
plt.figure;
plt.plot(mon, time_series)
plt.plot(mon, 0*time_series, 'k')
plt.xlim(init_year, fin_year)
plt.ylim(lower, upper)
plt.show()
フーリエ解析とは¶
前回は「フーリエ級数」を勉強しました。一定の周期を持つ時系列を，三角関数の足し合わせで表現する方法でした。今回勉強する （離散）フーリエ解析 では，気象や気候のデータ解析はもちろん電気信号処理や音楽情報処理などで頻出する「周期を持たない，有限の，離散的な時系列」にも使えるように拡張します。
ただ，この拡張の仕方をきちんと数学的に扱おうとすると，最低でも半期分くらいの授業を使うことになります。そこで今回は，数学的なことやアルゴリズム*に興味のある方には自分で勉強してもらうことにして，実際にPythonの組込関数を使って気候データに当てはめてパワースペクトルを計算してみることで，先行研究のデータを読んだり，ご自身の研究に当てはめるためのきっかけを得ること，くらいを目標にします。
- 特にアルゴリズムについては，本当は情報科学出身の方には全員知って欲しいくらい美しい高速フーリエ変換 Fast Fourier Transform (FFT) というアルゴリズムを用います。名前だけでも覚えて帰ってください。
Niño3.4指数のパワースペクトル¶
Niño3.4というインデックスがあったのを覚えていますでしょうか。各時刻の熱帯太平洋海面水温が，エルニーニョ寄りなのかラニーニャ寄りなのかを表す時系列です。
nino34 = aave(190, 240, -5, 5)
plot_mon_time(nino34)
この時系列を見て，おおよそどのくらいの時間スケールでエルニーニョとラニーニャを行ったり来たりしているのかを知りたいというのが，フーリエ解析のモチベーションです。
そこで，この時系列全体を1周期として繰り返しているような周期関数を考えて，そのパワースペクトルを求めるということをします。
（そのため，始点と終点が離れすぎていると，周期のつなぎ目の影響が結果に大きく反映されてしまい，周期関数への拡張があまり具合の良くないものになってしまいます。たとえば明確な上昇・下降トレンドがある場合は除去すべきです。）
プログラムとしては，実は前回のようにわざわざ三角関数の回帰係数をfor文で計算しなくても，numpyモジュールには高速フーリエ変換(FFT)の組み込み関数が与えられているので，以下のようにすれば良いです（まずは写経でOK，少しずつ勉強していってください）。
※神山はこちらも参考にしました。https://jp.mathworks.com/help/matlab/math/basic-spectral-analysis.html
N = nino34.shape[0] # サンプルサイズ
delta_t = 1 # 月別データのサンプリング間隔は1ヶ月
frequency = np.fft.fftfreq(N, delta_t)[0:math.floor(N/2)] # 計算する各周波数の値が入った配列
nino34_k = np.fft.fft(nino34)[0:math.floor(N/2)] # sinとcosの振幅のようなもの（D問題）
power = np.abs(nino34_k)**2*delta_t/N # パワースペクトルは振幅の二乗に比例
# 発展的な注：delta_t/Nをかけるのは，パワースペクトルを積分した時に全分散と一致させるため
##（厳密には，ここでは「パワースペクトル密度」を求めている）
plt.plot(frequency, power, 'r')
plt.xlabel("Frequency (/month)")
plt.ylabel("Power ($°C^2 \mathrm{month}$)")
plt.show()
ここで，横軸の刻み幅はサンプリング間隔を$\Delta t$，サンプルサイズを$N$とすると，$\frac{1}{N \Delta t}$となります。
発展的な注：このデータだと38年周期の逆数ですので，ちょうどデータの長さが38年分であることと対応しています。これよりも長い周期について知りたい場合は，元のデータを長くするしかありません。
実際上は，これほど細かい刻み幅でスペクトルを知らなければいけないことは稀です。
そこで，隣接するいくつかの周波数同士を平均することによって統計的自由度をかせぎ，ノイズを減らすことが大事です。
（あるいは，時系列をいくつかの短い時系列に分割して計算したそれぞれのパワースペクトルの平均を求めることによってノイズを減らす方法もあります）
# ここでは，隣接する4つの周波数のパワーをそれぞれ平均して「平均のパワー」を出すことにする。
# 自分で好きに決めて良い数字だが，
# あまりに大きすぎると横軸が大雑把すぎてどの周波数に注目すればいいかわからないし，
# あまり小さすぎるとノイズが減らないので，うまい塩梅を見つけるのが重要。
ave_num = 4
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]);
plt.plot(frequency_mean, power_mean, 'r')
plt.xlabel("Frequency (/month)")
plt.ylabel("Power ($°C^2 \mathrm{month}$)")
plt.show()
いまは理論的に求められるすべての周波数帯について描画していますが，もう少し注目したい周波数帯だけ描画することにしましょう。
plt.plot(frequency_mean, power_mean, 'r')
plt.xlim(0, 0.1)
plt.xlabel("Frequency (/month)")
plt.ylabel("Power ($°C^2 \mathrm{month}$)")
plt.show()
こうすると，0.02/month くらいの周波数帯（逆数をとると約50ヶ月周期）にピークがあることがわかります。
つまり，エルニーニョ南方振動現象は，4年程度の周期に最も大きなパワーを持った振動現象（＝4年に1度の周期でエルニーニョとラニーニャを繰り返すくらいの振動現象）であることがわかるのです。
東京の気温のパワースペクトル¶
さっそく他にも色々なデータで，パワースペクトルを計算してみましょう。統計的性質が一定であるような時系列ほどうまくいきます。
その前に，時系列を入れるとパワースペクトルを描いてくれるような関数を作っておきます。
# 引数は，時系列データ，サンプリング周期，いくつの周波数を平均するかの数
def draw_power(time_series, delta_t = 1, ave_num = 4):
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*delta_t/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]);
plt.plot(frequency_mean, power_mean, 'r')
plt.xlabel("Frequency (/month)") # 月別データ以外の場合は単位を書き換える
plt.ylabel("Power ($°C^2 \mathrm{month}$)") # 温度のデータ以外の場合は単位を書き換える
plt.show()
こうしておくことで，Niño3.4指数のパワースペクトルは1行で計算できるようになります。
draw_power(nino34)
それでは，東京の気温の時系列と気温偏差の時系列ではどうでしょうか。先ほど注意した通り，事前にトレンドを除去しておきましょう。
東京の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv （出典は環境情報論第3回に記述）
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]
# デトレンド（フーリエ解析では絶対に必要！！！）
tokyo = signal.detrend(tokyo)
tokyoa = signal.detrend(tokyoa)
それでは，早速東京の気温について，時系列とパワースペクトルを描いてみましょう。
plot_mon_time(tokyo, lower = -15, upper = 15, init_year=1990, fin_year=2019)
この時系列を見て，どんなパワースペクトルになるのか想像できますでしょうか？
draw_power(tokyo, 1, 4)
0.1/month のちょっとだけ左に，ものすごく大きなピークが出ていますね。
これは，季節変動（1/12 /month，つまり1年周期）がとても大きなデータであることの表れです。
次に，東京の気温偏差について，時系列とパワースペクトルを描いてみましょう。
plot_mon_time(tokyoa, lower = -4, upper = 4, init_year=1990, fin_year=2019)
draw_power(tokyoa, 1, 10)
気候値を引き算するだけでは完全に季節変動が除去できていないのか，少しだけ一年周期の変動が残っていますが，それでもだいぶなめらかなパワースペクトルになりました。
全体的に，低周波側に大きなパワーがあるのはなぜでしょうか？→用語解説へ
用語解説：レッドノイズとホワイトノイズ¶
実は，パワースペクトルを左端から右端まで積分してグラフの下の面積を求めると，元の時系列の分散に一致します。それゆえ，パワースペクトルの下のある一部分の面積に注目すれば，特定の周波数帯が全体に対してどの程度のパワーを持つか（元の分散に対してどの程度の寄与率があるか）を知ることができます。
たとえば，先ほどから見てきたような海面水温偏差や気温偏差のデータは，低周波（長周期）成分側にパワーが偏っています。 そのような性質をもつ時系列をレッドノイズといいます。これに対して，低周波から高周波までパワーの偏りがないような時系列をホワイトノイズといいます。
ある時刻のデータが，それよりも前の時刻のデータの影響を受けて決まる場合はレッドノイズになり，各時刻のデータが独立に決まる場合はホワイトノイズになります。たとえば，先月が暖かければ今月も暖かい可能性はそれなりに高いですから，気温はレッドノイズになります。海面水温だと，過去の水温と今月の水温の関係は気温よりもさらに強いですから，さらにレッドノイズ的になります。
気象は多くの場合何らかの形で過去を引きずるので，気象データは多かれ少なかれレッドノイズですが，たとえば雨の音を録音した波形の時系列は典型的なホワイトノイズです（面白いので色々と例をググってみてください）。
低周波（長周期）成分側にパワーが偏っている（＝強いレッドノイズである）ほど，その系は過去の情報を覚えているとか，過去のメモリーを持っているといいます。
気候データの中では，海洋から得られるデータの方が，大気のデータよりも概して強いレッドノイズになりがちです。実は，大気のメモリーは多くても数週間から数ヶ月ですが，海洋のメモリーは数年程度以上あります。 それゆえ，長期予報を行うためには海洋内部の情報が欠かせません。
また，海氷面積の時系列も強いレッドノイズですので，海氷も過去の情報をかなり覚えていると言えます。真冬で辺り一面に張っていた氷がある日突然目の前から消えたり，真夏で一週間前まで全く氷が見えなかった状態からある日急激に氷に覆い尽くされることはないということです。
課題2（締切：次回の授業開始前まで）¶
この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。
A問題. パワースペクトルとは何でしょうか？自分の言葉で説明してください。
B問題. draw_power関数の引数であるave_numを変えると，描画されるパワースペクトルがどのように変わるのかをいくつか実験し，簡単に考察してください。
C問題. あなたの趣味に関するデータなど，好きなデータのデトレンドした偏差についてパワースペクトルを計算し，レッドノイズ的か，ホワイトノイズ的かを考えてみましょう。
D問題. 次の二つについて答えてください。
I. フーリエ級数 $$ \begin{eqnarray} f(x) &=& \frac{a_0}{2} + \left\{a_1 \cos(\frac{\pi x}{L}) + b_1 \sin(\frac{\pi x}{L})\right \} + \left\{a_2 \cos(\frac{2 \pi x}{L}) + b_2 \sin(\frac{2 \pi x}{L})\right \} + ... \\ &=& \frac{a_0}{2} + \sum_{n=1}^\infty \left\{a_n \cos(\frac{n \pi x}{L}) + b_n \sin(\frac{n \pi x}{L})\right \} \end{eqnarray}$$ にオイラーの公式を用いると，三角関数の代わりに指数関数を用いたフーリエ級数（「複素フーリエ級数」）の公式 $$ \begin{eqnarray} f(x) &=& \sum_{n=-\infty}^\infty c_n e^{i \frac{n \pi x}{L}} \end{eqnarray}$$ が導出できます。$c_n$と$c_{-n}$（ただし$n = 0, 1, 2, ...$）を$a_n$, $b_n$を用いて表すとそれぞれどう書けるでしょうか？
II. nino34_k = np.fft.fft(nino34)では，各$n$についてNiño3.4の$c_n N$を求めています（$N$はサンプルサイズ）。そして，その次の行のpower = np.abs(nino34_k)**2*delta_t/Nでは，$|c_n|^2N\Delta t$を求めています。$|c_n|^2$を$a_n$, $b_n$を用いて書き直すことで，このようにして求めたパワースペクトルが，前回の授業で説明したように「それぞれの周波数についてのsinとcosの振幅の2乗（に比例する量）」を表していることを確かめてください。
発展的な注：パワースペクトルの横軸の刻み幅は$1/N\Delta t$ゆえ，パワースペクトルを全周波数について積分すると $$\sum_{n = 0}^{N-1} (|c_n|^2N\Delta t) (1/N\Delta t) = \sum_{n = 0}^{N-1} |c_n|^2$$ となる。パーセバルの定理より，これは元の時系列の分散 $$\frac{1}{N}\sum_{t=0}^{N-1} \{f(x_t)\}^2$$ に一致する。すなわち，パワースペクトル（密度）の下の面積は，その範囲の周波数帯が説明する分散の寄与率である。
提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論2，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。