ある数量の時系列が、ある周期τ、たとえば1年でもとにもどる周期現象であれば、 平均値と、 周期1年、1/2年、1/3年...の三角関数(sin, cos)の重ね合わせ(線形結合)で 表現することができる。式でかけば
∞
x(t) = a(0) + Σ { a(k) cos(2πkt/τ) + b(k) sin(2πkt/τ) } (1)
k=1
のようになる。
春日(1993)の式(4.2)は、 (1)式のa(0)のところを、a0/2と書いている。 つまり、春日のa0はここでのa(0)の2倍である。 「実数のフーリエ級数」の定義としては、春日の(4.2)の形が標準的であるが、 ここでは、a(k) cos(2πkt/τ) で k = 0にしたとき( cos(0) = 1 )に対応すると 考えて、(1)の形を採用しておく。
観測データでは、時系列に関する情報は、時間の連続関数で得られているわけ ではなく、多くの場合、一定の時間間隔でとられた離散的な時点での値が 得られている。 1年周期現象であれば、月ごとの値をもとに論ずることが多い。 現実の1か月の長さは一定でないが、ここでは等間隔とみなして話を進める。 有限個の観測値から、(1)式の無限個の係数の情報を得ることは不可能である (形式的には計算できるが、相互に独立でなくなる)。 n個の時点 tj = (j-1)Δt (j = 1,...,n) での値xjを表現するには、次のように、平均値のほかに (n/2)個のcosの係数と(n/2-1)個のsinの係数があればじゅうぶんである。
n/2 n/2-1
xj = a(0) + Σ a(k) cos(2πktj/τ) + Σ b(k) sin(2πktj/τ) (2)
k=1 k=1
1年周期現象を月値で表現する場合ならば、年平均値のほか、
(1/6)年周期までのcos成分と(1/5)年周期までのsin成分があればよい。
ここで、(1/6)年周期は、このサンプル間隔で表現できるかどうかの
ぎりぎりの振動数(信号処理の専門ではNyquist振動数という)に対応する。
波形を具体的に考えてみると、t = 0を含む等間隔のサンプリングでは、
(1/6)年周期のsin成分は常に0になってしまい、係数が決められない
ことがわかる。
{ xj } から { a(k), b(k) } を求めるような変換を、 離散(discrete)フーリエ変換ということがある。 なお、高速フーリエ変換(略称FFT)は、 この離散フーリエ変換を効率よく行なうアルゴリズムをさす。
この係数を求める方法はいろいろあるが、 「対象となる現象に周期性があり、 その周期の整数分の1である時間間隔でデータが得られており、 欠けがない」という条件が満たされているときは、 基底関数系となる三角関数の間の 「直交性」(違うものをかけて1周期積分すると0になること)を利用し、 もとの時系列にそれぞれの三角関数をかけてやることによって求めるのが いちばん簡単である。 時系列{xj, j = 1,...,n}が(2)のように表現されたと仮定し、 両辺に(1)の各基底関数をかけて、 1周期(周期で無次元化した時刻 (j-1)/n が 0 から 1 まで)積分してやればよい。 基底関数の違うものどうしの積の積分が消えるので、結果は簡単な形になる。 基底関数の同じものどうしをかけて1周期積分して得られる値は、 1とは限らないので、係数を求める場合はこれに注意して適切に重みづけ (「規格化」)する必要がある。
以上の計算を、n=12の場合について、 実際にためしてみる例文(Fortran)が harmonic.f, である。 標準入力から1行に1つずつ12個の数値を読みこみ、 まず、12個の係数を計算して(サブルーチンfourier)、書き出す。 それから、この係数にもとづいて、 12個の等間隔の時点での時系列の値を再合成し (サブルーチンsynthesis)、 入力された値とならべて書き出す。 係数を求める部分の計算は、次のようになっている。 ただし、便宜上、入力データの時間間隔を1とし、 基本周期τがn(=12)に等しいとしている。
1 n
a(0) = - Σ xj
n j=1
は時間平均そのものである。
2 n
a(k) = - Σ xj cos(2π(j-1)k/n)
n j=1
2 n
b(k) = - Σ xj sin(2π(j-1)k/n)
n j=1
1 n
a(k) = - Σ xj cos(2π(j-1)k/n)
n j=1
となる。
n/2 n/2-1
xj = a(0) + Σ a(k) cos(2πk(j-1)/n) + Σ b(k) sin(2πk(j-1)/n)
k=1 k=1
となる。
同じ計算をC言語で書いた例文が harmonic.c である。 ただし、次のような細かい違いがある。
Unix上でこのプログラムを動かしてみる手続きの例文が tharmonic_sh である。 入力データ例 kion.tokyo.1979.txtは、 周期現象の例としてはあまりよくないが、 特定の年(1979年)の東京の月平均気温である。 このデータファイルは、各行に、空白区切りで1列めに時間座標(年の小数値)、 2列めに気温がはいっているので、 awkで2列めを取り出して、harmonic.cをコンパイルした プログラムの標準入力に与えている。
ここまでは、時間の原点 t = 0 を、最初のデータ値が与えられている時点に とっていた。 しかし、データが月平均の量(たとえば月平均気温)の場合、最初のデータ値は 1月の気温だが、時間の原点は、1月のまんなかではなく、1月の初めにとるべき だろう。 そうすると、これまでに求めた三角関数成分は、すべて位相が0.5/12周期、 つまりπ/12ラジアンあるいは15°だけずれることになる。 この点を修正した例文(Fortran)が harmoni2.f である。 (なお、(1/6)年周期成分にcos, sinの両方が出てくるが、 係数は相互に独立ではない。)
このようにして求められた係数が、 1年の気温変化をなめらかにしたものを表わすと考えて、 細かい時間ステップ(たとえば毎日)の気温を求めることができるが、 これを月ごとに平均すると、入力として使った月平均気温と合わない。 実際に、毎日の気温を求めて月平均を計算してみるプログラム(Fortran)が harmoni3.f である。ただし、各月の日数をすべて30日とした。
月平均が合うような毎日の気温を求めるための係数は、 月平均気温を三角関数に分解して得られた係数に 補正を加えることによって得ることができる。 坪井(1968)のp. 254 - 256で、その補正の考え方が説明されている。 その結果をとり入れたプログラム(Fortran)が harmoni4.f であり、実際数値を入れてみると、ほぼ入力に一致した月平均値が得られる。 多少の誤差は、時間を連続で計算した場合に正しい補正をしながら、 実際には1日きざみで計算していることからくる。