数値解析法及び演習第十四回

[1] フーリエ変換の復習

フーリエ変換はさまざまな場面で出会うことになる. 地震の波形解析, 気候変動の解析, 地層の縞の解析などが典型的な例である. 今週と来週は離散データをフーリエ変換する手法を学ぶ. 今週は離散フーリエ変換のアルゴリズムを作成し, 来週はそれを用いる際の注意点および応用を学ぶ.

時間の関数h(t)が与えられたとしよう. h(t)のフーリエ変換, その逆変換は次式で与えられる. （注：一般には2πfではなくてωが用いられるが, ここでは計算の便利のため2πfを使う）

$\begin{displaymath} H(f)=\int^{\infty}_{-\infty}h(t)e^{2\pi ift}dt \end{displaymath}$

(1)

$\begin{displaymath} h(t)=\int^{\infty}_{-\infty}H(f)e^{-2\pi ift}dt \end{displaymath}$

(2)

フーリエ変換を行うと関数h(t)は関数H(f)に変換される. H(f)は一般に複素数となる. その絶対値の二乗は, 周波数f～f+dfに含まれるパワーになる. したがって, フーリエ変換すると, どの周波数成分が卓越しているかが分かる事になる. また，h(t)が実数の場合はH(-f)=H(f)^*(H(f)の複素共役）となることに注意．

以下でフーリエ変換を議論するときは, 次の形式でデルタ関数を書いておくと都合がよい.

$\begin{displaymath} \delta (t)=\int_{-\infty}^{\infty}e^{2\pi i ft}dt \end{displaymath}$

(3)

練習問題1

式（3）を用い, 式（1）から逆変換の式（2）を導け.

[2]離散フーリエ変換

我々が得るデータは普通離散データである. この場合, 一定の時間間隔Δtで物理量が計測される. 次の時刻t_nについてN個のデータが得られているとしよう.

$\begin{displaymath} t_n={n\Delta t},\quad 0\le n\le N-1 \end{displaymath}$

(4)

Δtについて, ナイキストの臨界周波数という次式で定義される周波数が存在する.

$\begin{displaymath} f_{\rm c}={1\over 2\Delta t} \end{displaymath}$

(5)

ナイキストの臨界周波数がなぜ重要かというと, この周波数以上の成分は正しくフーリエ変換されないためである. このため, 離散フーリエ変換をする時の周波数の上限が決まる.

ナイキストの臨界周波数の存在から, 離散フーリエ変換を行った場合に我々が得ることのできるフーリエ変換成分は-fcからfcまでであることが分かった. データがN個あるものとすると, このデータから得られるフーリエ変換成分もたかだかN個しかない. そこで次の値についてフーリエ変換を求めることにしよう.

$\begin{displaymath} f_n={n\over N\Delta t},\quad -{N\over 2}\le n\le {N\over 2} \end{displaymath}$

(6)

式（1）を和で置き換えると以下のようになる.

$\displaystyle H(f_n)$	$\textstyle =$	$\displaystyle \int^{\infty}_{-\infty}h(t)e^{2\pi if_nt}dt\simeq \sum_{k=0}^{N-1}h_ke^{2\pi if_nk\Delta t}\Delta t$
	$\textstyle =$	$\displaystyle \Delta t\sum_{k=0}^{N-1}h_ke^{2\pi i kn/N}, \quad -{N\over 2}\le n\le {N\over 2}$	(7)

一般に離散フーリエ変換は上式からΔtをとった

$\begin{displaymath} H_n=\sum_{k=0}^{N-1}h_ke^{2\pi i kn/N} \end{displaymath}$

(8)

を指す．このように定義しておくと次元によらずまた逆変換との対応もよく見える．

練習問題2

H_-n=H_N-nが成立している事を確認せよ. ここから，Hは周期Nの周期関数であることになる．

上の問題から, H_-n=H_N-nであることが分かる. もともと式（4 ）においてnは-N/2からN/2まで変化していたが, 0からN-1まで変化させることにしよう. こうすると, 0からN/2までは周波数0からfcまで対応し, N/2からN-1までは周波数-fcから-1/NΔtまでが対応することになる. すると最終的に求めるべき離散フーリエ変換は次式となる.

$\begin{displaymath} H_n=\sum_{k=0}^{N-1}h_ke^{2\pi i kn/N}, \quad 0\le n\le N-1 \end{displaymath}$

(9)

n番めの周波数のパワーは絶対値の二乗

$\begin{displaymath} \vert H_n\vert^2 \end{displaymath}$

(10)

で与えられる．もとの関数が実数である場合には，周波数fと周波数-fのパワーは一致することに注意．

また, 離散逆フーリエ変換は次式で与えられる.

$\begin{displaymath} h_k={1\over N}\sum_{n=0}^{N-1}H_ne^{-2\pi i kn/N}, \quad 0\le n\le N-1 \end{displaymath}$

(11)

練習問題3

$\displaystyle \sum_{k=0}^{N-1}e^{2\pi i k(n-m)/N}$	$\textstyle =$	$\displaystyle N \quad n=m$
	$\textstyle =$	$\displaystyle 0 \quad {\rm else}$	(12)

を用いて, 式（9）から離散逆フーリエ変換の式（11）を導け.

１）式（9）を元に, 離散フーリエ変換を行うアルゴリズムを作成せよ. もとのデータは実数とする.

２）式（11）を元に, 離散フーリエ変換の結果を逆変換するアルゴリズムを作成せよ.

[3]様々な関数の離散フーリエ変換

このセクションでは，作成したプログラムを用いて，様々な関数の離散フーリエ変換を行う．

定数

このデータは，時間間隔 0.209439(=2π/30)，時間領域 0< t < 2π において定数関数h(t)=1のデータである．データ数は３０である．このデータをセーブしてgnupotでグラフを作成せよ．離散フーリエ変換を行い

どの周波数のパワーが強くなっているか

フーリエ成分の実部と虚部の大きさをgnuplot でプロットすることによりそれぞれ確認せよ．式（９）と（１２）からhk=1の場合の離散フーリエ変換を行ってみよ．数値計算の結果と比較せよ．

cos関数とsin関数

h(t)=cos(t)のデータは，時間間隔 0.209439(=2π/30)，時間領域 0< t < 2π においてcos(t)の値を出力したものである．データ数は３０である．このデータをセーブしてgnupotでグラフを作成せよ．離散フーリエ変換を行い

どの周波数のパワーが強くなっているか

フーリエ成分の実部と虚部の大きさをそれぞれ確認せよ．また，h(t)=sin(t)のデータについても行うこと．

sin関数その2

次に，sin(t)から位相がずれてsin(t+Δ)となっている場合の離散フーリエ変換をおこなってみよ． h(t)=sin(t+Δ)のデータは，先ほどと同じ時間間隔，時間範囲においてsin(t+Δ)の値を出力したものである．

データをgnuplotで確認し
離散フーリエ変換を行い，位相Δを決定せよ．

sin関数その3

このデータは，3つの三角関数の和を先ほどと同じ時間間隔，時間範囲で出力したものである．離散フーリエ変換を行い，3つの周波数を決定せよ．

データ１は，h(t)=sin(20t)を時間間隔 0.209439(=2π/30)でサンプリングしたものである．この場合のナイキストの臨界周波数は，fc=1/2Δt=15/2πとなる．sin(20t)の周波数は20/2πであるので臨界周波数よりも大きい．このデータの離散フーリエ変換を行い，どの周波数が強くなっているか確認せよ．またデータ２は，h(t)=sin(50t)をサンプリングしたものである．gnuplotでデータ１と重ねてプロットしてみよ．

有限の時間範囲

離散フーリエ変換は当然有限の時間範囲のデータを取り扱うことになる．ここから引き起こされる現象を確認する．次の二つのデータは，時間間隔 Δ t=0.209439(=2π/30)で0&le t&le 999*Δtの範囲の1000個のサンプルである．一方は初めの10個が1，残りすべては0であり，もう一方は初めの100個が1で残りはゼロである．

二つのデータをgnuplotで確認せよ．
離散フーリエ変換を行い，パワーの分布をgnuplotで確認せよ．

先ほどと同じ時間間隔，時間範囲でh(t)=1/tのグラフを100*Δtまでサンプルしたものと10*Δtまでサンプルしたものをそれぞれ離散フーリエ変換し，その結果を比較せよ．この結果から，サンプリングする時間範囲を長くすればするほどリップルの効果が小さくなることが分かる．

今日の宿題 (1/27日まで)

実際に得られたデータのフーリエ変換を行ってみる．データに含まれている強い波の周期を読みとってみよ．

南極ボストーク基地の二酸化炭素濃度（平均値からのずれ）：約42万年間, 500点, Δt=825.136年

プログラムと上のデータに含まれる強い周期（年）を送ること. その周期は地球の軌道変化の周期と関連していると言われている．

日程表へ戻る << 　

数値解析法及び演習 第十四回