離散 Fourier 変換
桂田 祐史
2014 年 11 月 14 日, 2014 年 11 月 21 日
応用上重要な離散 Fourier 変換の定義と簡単な性質を紹介する。用いる数学は複素数と線形
代数である (有限個のデータの話で、収束の話 (解析学) は必要ない)。むしろ Fourier 変換よ
りも簡単と言える。
離散 Fourier 変換には、高速なアルゴリズムである高速 Fourier 変換 (Fast Fourier Transform,
FFT) があるが、詳細は省略する。
概要
離散 Fourier 変換は、Fourier 係数、Fourier 級数展開の離散化と言える。
周期 2π の周期関数 f があるとき、測定して得られるような離散的な標本データを用
いて、近似的に Fourier 級数展開をすることを考える。適当な自然数 N を用いて、周期区
∫ π
1
間を N 等分し、N 等分点上の関数値 {fj } から Fourier 係数 cn =
f (x)e−inx dx
2π −π
の近似値 Cn を求めることになるが、それには積分を台形則で近似するのが良い。そうし
て得られた {Cn } は周期 N の周期数列である。写像
f0
C0
f1
C1
N
C ∋ . 7→ . ∈ CN
..
..
fN −1
CN −1
1 ( −(n−1)(j−1) )
ω
(ただし ω = e2πi/N )
は、離散 Fourier 変換と呼ばれるが、それは W :=
N)
√
(
を表現行列とする線型写像で、W −1 = ω (j−1)(n−1) である。 N W は unitary 行列に
なっている (すなわち離散 Fourier 変換は実質的に unitary 変換である)。
1
離散 Fourier 係数 — なぜそのように定義するか
(以下、周期を 2π として説明するが、π を L, nx を nπx/L と書き換えることで、周期 2L
の場合の式が得られる。)
f : R → C が周期 2π で、ほどほどに良い関数とする (この文書の中で収束の問題は考えな
いことにする)。
∫ 2π
1
(1)
cn :=
f (x)e−inx dx (n ∈ Z)
2π 0
1
とおくとき (これまで [−π, π] で積分したが、[0, 2π] で積分しても同じ値が得られる)、
∑
(2)
f (x) =
cn einx .
n∈Z
1 周期区間 [0, 2π] を N 等分して、等分点
xj = jh (j = 0, 1, · · · , N ),
(3)
h :=
2π
N
での値
fj := f (xj )
(j = 0, 1, · · · , N )
を用いて 1 、{cn } の近似 {Cn } を求めよう。
周期関数の周期区間上の積分の数値計算には、台形公式を用いるのが良いとされる (やって
みると分かるが、しばしば驚異的な高精度が得られる)。
一般に定積分
∫
b
I=
F (t) dt
a
に対する台形公式とは、{tj }N
j=0 を [a, b] の N 等分点として、I を
TN :=
N
∑
F (tj−1 ) + F (tj )
j=1
2
(
′
′
h =h
F (t0 )
F (tN )
+ F (t1 ) + · · · + F (tN −1 ) +
2
2
)
,
h′ :=
b−a
N
で近似するものであるが、F が周期 b − a であれば、F (t0 ) = F (a) = F (b) = F (tN ) であるか
ら、次のような簡潔な式に帰着される:
′
TN = h
N
−1
∑
F (tj ) (周期関数に対する台形則).
j=0
(1) を台形則で近似したものを Cn とおこう。
1 ∑
h
f (xj )e−inxj .
2π j=0
N −1
Cn :=
(4)
ω := e2πi/N = eih
とおくと、
e−inxj = e−injh = ω −nj
であるから、
N −1
1 2π ∑
Cn =
fj ω −nj .
·
2π N j=0
j ∈ Z に対して xj , fj を考えることも出来るが、その場合、f が周期 2π であるので、{fj }n∈Z は fj+N = fj
(j ∈ Z) を満たす、すなわち周期 N の周期数列である。従って、{fj }j∈Z を扱う際は、連続した N 項、例えば
−1
{fj }N
j=0 だけを知れば十分である。つまり、j の範囲を広げることにあまり意味はない。
1
2
すなわち
N −1
1 ∑
Cn =
fj ω −nj .
N j=0
(5)
これを f の離散 Fourier 係数と呼ぶ。
補題 1.1 (1 の N 乗根の性質) N ∈ N に対して
ω = e2πi/N
とおくとき、次の (1), (2) が成り立つ。
(1) ω は 1 の原始 N 乗根、すなわち
(i) 1 ≤ m ≤ N − 1 ⇒ ω m ̸= 1
(ii) ω N = 1
を満たす。
(2) ∀m ∈ Z に対して、次式が成り立つ。
{
N
−1
∑
N (m ≡ 0 (mod N ))
mj
ω =
0 (それ以外).
j=0
証明
(1) k ∈ Z に対して
2π
ω k = ei N k
であるから、k ≡ 0 (mod N ) のとき、そのときのみ ω k = 1 が成り立つ。特に 1 ≤ k ≤ N −1
を満たす任意の k に対して、ω k ̸= 1, ω N = 1 であることに注意する。
(2) m ≡ 0 (mod N ) のとき、ω m = 1 であるから、任意の j に対して ω mj = 1. ゆえに
N
−1
N
−1
∑
∑
mj
ω =
1 = N . そうでないとき、ω m ̸= 1 であるから、等比数列の和の公式より
j=0
j=0
N
−1
∑
j=0
ω mj =
N
−1
∑
j=0
( )m
1 − ωN
1 − (ω m )N
(ω ) =
=
= 0.
1 − ωm
1 − ωm
m j
以下、この節では、2 数 ℓ, m がともに整数であり、ℓ ≡ m (mod N ) であることを単に ℓ ≡ m
と書くと約束する。
3
命題 1.2 (離散 Fourier 係数の性質) 周期 2π の周期関数 f : R → C と N ∈ N に対して、
h :=
2π
,
N
xj := jh,
ω := e2πi/N = eih ,
(j ∈ Z),
fj := f (xj )
N =1
1 ∑
Cn :=
fj ω −nj
N j=0
(n ∈ Z)
で離散 Fourier 係数 {Cn }n∈Z を定めるとき、以下の (1), (2) が成り立つ。
(1) {Cn }n∈Z は周期 N の周期数列である: Cn+N = Cn (n ∈ Z).
(2) 任意の n ∈ Z に対して、
(6)
∑
Cn =
cm .
m≡n
ただし
∑
は、m ≡ n (mod N ) を満たすすべての m に対して和を取ることを意味
m≡n
する。
証明
(1) ω −(n+N )j = ω −nj ω −N j = ω −nj であるから、
Cn+N
N −1
N −1
1 ∑
1 ∑
−(n+N )j
=
fj ω
=
fj ω −nj = cn .
N j=0
N j=0
(2)
fj =
∑
cn einxj =
∑
n∈Z
n∈Z
cn einjh =
∑
cn ω nj
n∈Z
であるから、
N −1
N −1
N −1
1 ∑
1 ∑
1 ∑ −nj
−nj
ω
=
fj ω
=
Cn =
N j=0
N j=0
N j=0
=
1 ∑
cm
N m∈Z
N
−1
∑
ω (m−n)j =
j=0
(
∑
)
cm ω mj
ω −nj
m∈Z
∑
1 ∑
cm N =
cm .
N m≡n
m≡n
−1
この命題の (1) から、{Cn }n∈Z を求めよ、と要求されたとき、連続した N 項、例えば {Cn }N
n=0
を計算出来れば十分である (それ以上やることはない) ことが分かる。
4
C0 =
∑
cm = c0 + cN + c−N + c2N + c−2N + · · · ,
m≡0
C1 =
∑
cm = c1 + c1−N + c1+N + c1−2N + c1+2N + · · · ,
m≡1
C−1 =
∑
cm = c−1 + c−1+N + c−1−N + c−1+2N + c−1−2N + · · · ,
m≡−1
C2 =
∑
cm = c2 + c2−N + c2+N + c2−2N + c2+2N + · · · ,
m≡2
C−2 =
∑
cm = c−2 + c−2+N + c−2−N + c−2+2N + c−2−2N + · · ·
m≡−2
次の問を考えよう。
Q Cn は cn の良い近似になっているか？
(もともと cn の近似値を求めようとして、その定義式の積分を、台形則という数値積分公式で
近似したものが Cn であるわけだが、本当に Cn が cn に近いかは、確かめる必要がある。)
A ある意味で Yes である。正確には、
任意の n ∈ Z を固定するとき、
lim Cn = cn .
N →∞
注意: h, xj , fj , ω, Cn は、N に依存していることに注意せよ。本来は、それぞれ hN , xj,N ,
fj,N , ωN , Cn,N と書くべきものかもしれない。だから上の式は lim Cn,N = cn である。すな
N →∞
わち、ε-N で表せば
(∀n ∈ N)(∀ε > 0)(∃n′ ∈ N)(∀N ∈ N : N ≥ n′ )
|Cn,N − cn | < ε.
次に考えてもらいたいのは、
Q C−1,N = CN −1,N = c−1 + c−1+N + c−1−N + c−1+2N + c−1−2N + · · · であるわけだが、これ
は c−1 , cN −1 , いずれの近似か？
A c−1 の近似になっているが、cN −1 の近似ではない。
lim (c−1+N + c−1−N + c−1+2N + c−1−2N + · · · ) = 0
N →∞
であるから、確かに lim C−1,N = c−1 .
N →∞
結局、Cn (|n| < N/2) は cn の近似である、とまとめられるであろう。|n| ≒ N/2 のときは、
精度はかなり悪くなる。
2
離散 Fourier 変換
ここからは N 次元ベクトルの変換の話で、無限和は登場せず、線形代数の話になる。
5
f =
C=
f0
f1
..
.
N −1
1 ∑ −nj
N
ω fj で定義される
∈ C に対して、離散 Fourier 係数の式 Cn =
N j=0
fN −1
C0
C1
∈ CN を、f の離散 Fourier 変換と呼ぶ。また写像 CN ∋ f 7→ C ∈ CN も離
..
.
CN −1
散 Fourier 変換と呼ぶ。これは線型写像であり、ある行列をかけることで表現できる。その
行列の性質を述べよう。
命題 2.1 (離散 Fourier 変換の表現行列とその逆行列) N ∈ N に対して ω := e2πi/N ,
ω0
ω0
ω0
···
ω0
0
ω −1
ω −2
···
ω −(N −1)
ω
1
1 ( −(n−1)(j−1) )
−2
−4
−2(N −1)
ω 0
ω
ω
·
·
·
ω
=
ω
W :=
,
N
N
..
..
..
..
...
.
.
.
.
ω 0 ω −(N −1) ω −(N −1)2 · · ·
とおくとき、
f =
f0
f1
..
.
C=
,
ω −(N −1)(N −1)
C0
C1
..
.
CN −1
fN −1
に対して
(7)
Cn =
N −1
1 ∑
fj ω −nj
N j=0
(n = 0, 1, · · · , N − 1)
⇔
C = Wf.
W は正則で、W −1 の (j, n) 成分は ω (j−1)(n−1) である。言い換えると、
ω0 ω0
ω0
···
ω0
0
ω1
ω2
···
ω N −1
ω
0
2
4
2(N −1)
ω
ω
ω
·
·
·
ω
W −1 =
.
..
..
..
..
...
.
.
.
.
0
N −1
(N −1)2
(N −1)(N −1)
ω ω
ω
··· ω
すなわち
Cn =
N −1
1 ∑
fj ω −nj
N j=0
(n = 0, 1, · · · , N − 1)
⇔
fj =
N
−1
∑
ω jn Cn
(j = 0, . . . , N − 1).
n=0
6
証明 前半の (7) は行列の積の定義からすぐ分かる。
以下、この証明の最後まで、行列の添字は 0 から始める (N − 1 までということになる) と
いうルールで記述する。
1
W の (n, j) 成分は ω −nj である。
N
この W に、(j, k) 成分が ω jk である行列 (ω jk ) を右からかけた行列の (n, k) 成分は
{
{
N
−1
N −1
∑
1 ∑ (k−n)j
1 −nj jk
1
N (k − n ≡ 0)
1 (k = n)
=
= δkn .
ω ω =
ω
=
N
N k=0
N
0 (それ以外)
0 (それ以外)
k=0
( )
これは積が単位行列であることを示す。すなわち W −1 = ω jk .
余談 2.1 既に fj =
N
−1
∑
cn ω nj であることは示してあるので、絶対収束を仮定すれば、和の順
n=0
序を交換して
fj =
N
−1
∑
∑
cm ω
mj
=
n=0 m≡n
N
−1
∑
∑
cm ω
nj
n=0 m≡n
=
N
−1
∑
ω
nj
n=0
∑
cm =
m≡n
N
−1
∑
ω nj Cn
n=0
と出来て、これから W −1 の成分 ω jn であると分かるが、W の性質は関数 f とは関係ないの
で、上の証明の方が良いであろう。
注意 2.2 (ちょっと修正すると unitary 行列になる) U :=
(8)
1 ( −(n−1)(j−1) )
,
U=√
ω
N
√
N W とおくと、
1
1 ( (j−1)(n−1) )
U −1 = √ W −1 = √
ω
N
N
となる。ω = ω −1 に注意すると、U の Hermite 共役 U ∗ は
1 ( (j−1)(n−1) )
1 ( −(j−1)(n−1) )
∗
= U −1 .
ω
=√
ω
U =√
N
N
すなわち U は unitary 行列である。
通常は Fourier 級数展開は
∫ 2π
1
cn =
f (x)e−inx dx,
2π 0
f (x) =
∑
cn einx
n∈Z
と定義されるが、
1
cn = √
2π
∫
f (x)e−inx dx,
0
{
と定義して (これは正規直交系
2π
√1 einx
N
1 ∑
f (x) = √
cn einx
2π n∈Z
}
で展開したことになる)、それに基づいて離散 Fourier
係数を定めると、変換の行列として U が導かれる。
以上のように定義し直したい誘惑があるけれど、混乱する人が出て来ると思うのでやめてお
く。次のように命題としてまとめておくことにする。
7
系 2.3 命題 2.1 の行列を W とするとき、U :=
√
N W は対称な unitary 行列である。
余談 2.2 (離散 Fourier 変換の応用上の強力さを別にしても、個人的に面白いと感じるところ)
Fourier 級数展開
∫ π
∑
1
f (x)e−inx , f (x) =
cn =
cn einx
2π −π
n∈Z
と、Fourier 変換と反転公式 (Fourier 変換を逆 Fourier 変換すると元に戻る)
∫ ∞
∫ ∞
1
1
−ixξ
b
f (x)e
dx, f (x) =
f (ξ) = √
fb(ξ)eixξ dξ
2π −∞
2π −∞
の対応について既に話してあるが、離散 Fourier 変換とその反転公式 (離散 Fourier 変換を逆
離散 Fourier 変換すると元に戻る)
N −1
1 ∑
fj ω −nj ,
Cn =
N j=0
fj =
N
−1
∑
Cn ω nj
n=0
も良く対応している。むしろ、Fourier 級数よりもバランスの良い変換になっていて (f は周
期関数で {cn } は無限数列であるが、{fj } と {Cn } は周期数列になる)、むしろ Fourier 変換
に似ているとも言える。
実は Fourier 変換、逆 Fourier 変換は、L2 (R) における unitary 変換というものになってい
る (ますます似ていると感じる)。
ω n·0
ω n·1
ω n·2
..
.
φn :=
n·(N −1)
ω
とおくと、
f=
N
−1
∑
cn φ n
n=0
これは C
A
N
の直交系 {φn } による f の展開になっている。
Mathematica で離散 Fourier 変換
Mathematica の関数 Fourier[f ] は、f = (f0 , f1 , . . . , fN −1 ) に対して
N −1
1 ∑
√
fn ω nj
N n=0
(j = 0, 1, · · · , N − 1)
8
を計算する (これが Mathematica の離散 Fourier 変換の定義)。逆変換は InverseFourier[C ]
で行なう。
可能な場合は高速 Fourier 変換を使って計算されるので、効率が良くなる (効率をあげたい
ときは N の値に注意すべきである)。
この講義の離散 Fourier 変換の定義に一致する結果を得るには、次のようなオプションを指
定する。
Fourier[f,FourierParameters->{-1,-1}]
B
離散 Fourier 変換の流儀
伊理 [1], [2]
大浦 [3]
参考文献
い り まさお
[1] 伊理正夫：一般線形代数, 岩波書店 (2003), 伊理正夫, 線形代数 I, II, 岩波講座応用数学, 岩
波書店 (1993, 1994) の単行本化.
い り まさお
[2] 伊理正夫：線形代数汎論, 朝倉書店 (2009), 「一般線形代数」のリニューアル.
[3] 大浦：FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法, http://www.kurims.
kyoto-u.ac.jp/~ooura/fftman/ (1997∼).
9
© Copyright 2024 ExpyDoc
