離散フーリエ変換
数値計算に向けた解説書

シキノ ∗June 18, 2022
Contents
1 まとめ

3

2 連続の場合

4

2.1

2.2

フーリエ変換 (Fourier-Transform, FT) . . . . . . . . . . . . . . . . . . . .

4

2.1.1

本稿のフーリエ変換の定義 . . . . . . . . . . . . . . . . . . . . . . .

4

2.1.2

フーリエ変換に関係する簡単な式と特徴 . . . . . . . . . . . . . . .

5

2.1.3

フーリエ変換の例 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

畳み込み (Convolution) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2.1

畳み込みの定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2.2

畳み込みとフーリエ変換 . . . . . . . . . . . . . . . . . . . . . . . .

7

2.2.3

畳み込みの例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83 離散の場合
3.1

3.2
∗

10

離散フーリエ変換 (Discrete Fourier-Transform, DFT) . . . . . . . . . . . . 10

3.1.1

位置を無限区間から有限区間へ制限 → 波数の刻み幅 ∆k の決定 . . 10

3.1.2

波数を無限区間から有限区間へ制限 → 位置の刻み幅 ∆x の決定 . . 11

3.1.3

順方向/逆方向離散フーリエ変換の表記 . . . . . . . . . . . . . . . . 14

3.1.4

離散フーリエ変換の簡略化 . . . . . . . . . . . . . . . . . . . . . . . 16

3.1.5

本稿の離散フーリエ変換の定義 . . . . . . . . . . . . . . . . . . . . 21

3.1.6

離散フーリエ変換の例 . . . . . . . . . . . . . . . . . . . . . . . . . 23

離散畳み込み (Convolution) . . . . . . . . . . . . . . . . . . . . . . . . . . 29

https://slpr.sakura.ne.jp/qp/, https://twitter.com/sikinote

1

CONTENTS

CONTENTS

3.2.1

離散畳み込みの定義 . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.2.2

離散畳み込みと離散フーリエ変換 . . . . . . . . . . . . . . . . . . . 30

3.2.3

任意区間における離散畳み込みと離散フーリエ変換 . . . . . . . . . 31

3.2.4

離散畳み込みの例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

4 プログラム

34

4.1

4.2

Intel MKL を利用した離散 (高速) フーリエ変換 . . . . . . . . . . . . . . . 34
4.1.1

Intel MKL の離散フーリエ変換の定義 . . . . . . . . . . . . . . . . . 34

4.1.2

プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.1.3

MKL を利用した離散フーリエ変換の計算例 . . . . . . . . . . . . . 39

Intel MKL を利用した畳み込み . . . . . . . . . . . . . . . . . . . . . . . . 43
4.2.1

4.3

Intel MKL を利用した畳み込みの計算例 . . . . . . . . . . . . . . . 43

DFT の計算速度. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.3.1

DFT の高速化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.3.2

DFT の分点数に対する計算速度 . . . . . . . . . . . . . . . . . . . . 46

4.3.3

並列計算時の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

A フーリエ変換→逆変換で元に戻ることの証明

49

A.1 連続の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
A.2 離散の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
B DFT の具体例で考えた無理数倍について

50

C 任意区間の離散フーリエ変換プログラム

52

D 任意区間の離散フーリエ変換を利用する畳み込みプログラム55

2

1 まとめまとめ

1

フーリエ変換∫ ∞f (x)e−i2πkx dx
 f (k) =
∫−∞
∞



f (k)ei2πkx dk f (x) =離散フーリエ変換

N
−1
∑




f (kn ) =f (xm )e−i2πnm/N

(1.1a)

m=0
N
−1
∑



1



 f (xm ) = N(1.1b)−∞f (kn )ei2πnm/N(1.5a)(1.5b)

n=0

位置

位置

x = [0, X]
x = [−∞, ∞](1.2a)

xm = m∆x, (m = 0, 1, · · · , N − 1)

(1.6b)

X
1=N
K(1.6c)

k = [0, K],(1.7a)∆x =

波数(1.6a)

波数k = [−∞, ∞](1.3a)
kn = n∆k, (n = 0, 1, · · · , N − 1)

畳み込み(f ∗ g)(x) ≡

∫ ∞
−∞

K
1
= ,
N
X

(1.7c)

1, KX = N
N

(1.8)

∆k =
f (τ )g(x − τ )dτ

(1.4a)(f ∗ g)(k) ≡ F[(f ∗ g)(x)](k)= f (k)g(k)∆x∆k =(1.4b)(1.7b)f (xm ± X) = f (xm ), f (km ± K) = f (kn )
(1.9)
畳み込み

(f ∗ g)d (xn ) ≡∑f (xl )g(xn − xl ) (1.10a)l(f ∗ g)d (kn ) ≡ F[(f ∗ g)d (x)](kn )= f (kn )g(kn )3(1.10b)2 連続の場合

2

連続の場合

2.1

フーリエ変換 (Fourier-Transform, FT)

2.1.1

本稿のフーリエ変換の定義

本稿ではフーリエ変換を以下の演算と定義します 1 .

∫ ∞



f (x)e−i2πkx dx
 f (k) =
∫−∞
∞



f (k)ei2πkx dk f (x) =(2.1a)(2.1b)−∞

ここで, Eq. (2.1a) を順方向フーリエ変換, Eq. (2.1b) を逆方向フーリエ変換と呼びます.
多くの場合で x は位置（時間）, k は波数（周波数）と呼ばれます.
フーリエ変換は右辺の積分を省略して{

f (k) = F [f (x)] (k)(2.2a)

f (x) = F −1 [f (k)] (x)(2.2b)と演算子 F, F −1 を用いて表記したりします.
フーリエ変換は, 基底関数 φk (x) = ei2πkx への射影と考えることができます. つまりEq. (2.1b) とは,『適当な関数 f (x) があるとき, f (x) を φk (x) の f (k) 倍の和（積分）で書いてください.』という問題なわけです 2 . ここで, φk (x) にかかる係数は x に依存しない k の関数 f (k) であり, それは Eq. (2.1a) から求められる, という操作がフーリエ変換なわけです.

1

フーリエ変換の定義は分野によって異なります. 基本的には, その分野で何が目的を元に決められます.

数値計算の分野では計算量が少なくなることが素晴らしいことので, 積分記号の前に余計な係数が付かないこと等から, 指数関数の肩に 2π が付けられたものが使われます. 物理学（特に量子力学）では, 角周波数で表現すると式の見た目が良くなるので指数関数の肩に 2π が付かないものが使われます.
2
なかなか積分は和と同じということをイメージするのは慣れが必要なので, 人によっては離散フーリエ変換 (節 3.1) の方がイメージしやすいかもしれません.

4

2 連続の場合

2.1.2

2.1 フーリエ変換 (Fourier-Transform, FT)

フーリエ変換に関係する簡単な式と特徴• 逆方向フーリエ変換から順方向フーリエ変換を導く
逆方向フーリエ変換の結果だけから, 順方向フーリエ変換の結果を導くことができ′ます 3 . Eq. (2.1b) の左から e−i2πk x を掛けて, そのあと全位置空間 x で積分すれば,

f (k) について導くことができます. つまり,
∫ ∞
∫ ∞
∫ ∞
−i2πk′ x−i2πk′ xdxe
f (x) =dxe·
dkf (k)ei2πkx−∞∫ ∞ −∞∫−∞∞′
dxei2π(k−k )x
dkf (k)
=
−∞
∫−∞
∞
=
dkf (k)δ(k − k ′ )−∞′= f (k )(2.3a)(2.3b)(2.3c)(2.3d)となります.k ′ を改めて k と書けば,∫ ∞f (k) =

dxe−i2πkx f (x)(2.3e)−∞となり, 順方向フーリエ変換が導けます.

• フーリエ変換はユニタリー変換
フーリエ変換の重要な特徴の一つはユニタリー変換であることです. つまり, フー
リエ変換して逆変換すると元の関数に戻ります. 式で示せば,

F −1 Ff (x) = f (x),→F −1 F = 1ということです. 実際に考えてみますと
∫ ∞
′
−1
−1
F Ff (x) = F
dx′ f (x′ )e−i2πkx[∫ ∞]∫ ∞ −∞′′ −i2πkx′=
dk
dx f (x )eei2πkx
−∞
∫−∞
∫ ∞
∞
′
′
′
=
dx f (x )dk ′ ei2πk(x−x )
−∞
∫−∞
∞
=
dx′ f (x′ )δ(x − x′ )(2.4)

(2.5a)(2.5b)(2.5c)(2.5d)−∞= f (x)(2.5e)となります.
3

順方向から逆方向を導いてもどちらでも構いません

5

2 連続の場合

2.1.3

2.1 フーリエ変換 (Fourier-Transform, FT)

フーリエ変換の例

実際にフーリエ変換を Eq. (2.1) に従って実行してみましょう.

f (x) = cos(2π3x)(2.6)という, 波数 3 を持つ波をフーリエ変換してみます. フーリエ変換すると波数の成分が出てくるはずなので, 波数 3 が答えとして出てくることを期待します.
順方向フーリエ変換を行うと 4 ,∫ ∞dx cos(2π3x)e−i2πkx
f (k) =−∞
∫ ∞ei2π3x + e−i2π3x −i2πkx
=
dx
e
2
−∞
∫
]
1 ∞ [ i2π(3−k)x=dx e+ e−i2π(3+k)x
2 −∞1=[δ (k − 3) + δ (k + 3)]
2(2.7a)(2.7b)(2.7c)(2.7d)と求められます. つまり, cos(2π3x) をフーリエ変換したら k = ±3 に正の振幅を持ったピークが出てくる, ということです. フーリエ変換によって, 波形の波数に値が出てくる,ということが確認できました.

4

積分の文字 dx の順番を積分の直後に置いておりますが, 同じ意味です. 何個も積分記号が混在すると,どの変数がどの範囲か分からなくなるため, 積分記号の直後に置くことにしています.

62 連続の場合

2.2 畳み込み (Convolution)

畳み込み (Convolution)

2.2
2.2.1

畳み込みの定義

畳み込み（もしくは畳み込み積分）とは, 関数 f (x), g(x) が定義されているときに以下の積分(f ∗ g)(x) ≡

∫ ∞
∫−∞
∞

=
−∞

f (τ )g(x − τ )dτ

(2.8a)f (x − τ )g(τ )dτ

(2.8b)

を計算することです. 関数 f と関数 g の畳み込み, と言います.
畳み込みとは『関数と関数を混ぜ込んだ結果』みたいなイメージです. Eq. (2.8a) から

Eq. (2.8b) へは, x − τ を新たな変数として変数変換すれば導けます. τ の範囲は (−∞, ∞)が基本であり本稿はそれに倣いますが, 問題によっては有限範囲や 1 周期を考える場合も
あります.

2.2.2

畳み込みとフーリエ変換

畳み込み積分 (2.8a) をフーリエ変換された関数で表現してみましょう 5 .
条件と問題を整理すると以下の通りになります.
∫(f ∗ g)(x) = dτ f (τ )g(x − τ )∫ ∞f (x) =
dkf (k)ei2πkx∫ −∞∞g(x) =
dkg(k)ei2πkx(2.9a)(2.9b)(2.9c)−∞

Eq. (2.9b), (2.9c) を Eq. (2.9a) に代入すれば,∫ ∞(f ∗ g)(x) =dτ f (τ )g(x − τ )−∞[∫ ∞] [∫ ∞]∫i2πkτ
′′ i2πk′ (x−τ )=dτdkf (k)edk g(k )e
−∞
−∞
[
]
∫ ∞
∫ ∞
∫
′
′ i2πk′ xi2π(k−k′ )τ
=dk
dk f (k)g(k )e
dτ e
−∞
−∞
∫ ∞
∫ ∞
[
]
′ i2πk′ x′′
δ(k − k )=dk
dk f (k)g(k )e
−∞
−∞
∫ ∞
=
dkf (k)g(k)ei2πkx(2.10a)(2.10b)(2.10c)(2.10d)(2.10e)−∞5

なぜフーリエ変換しようと思ったかは聞かないでください. 先人が畳み込みを計算しようと思って考えついて, それが綺麗だった, 程度の理由だと思います.

7

2 連続の場合

2.2 畳み込み (Convolution)と簡単になります. つまり, ei2πkx の前にかかる係数がフーリエ変換結果ですので, 畳み込みの結果をフーリエ変換すると,

F[(f ∗ g)(x)](k) =∫ ∞−∞(f ∗ g)(x)e−i2πkx dx= f (k)g(k)(2.11a)(2.11b)となります.
以上より『畳み込みの計算結果をフーリエ変換したい！』となったら, 畳み込みに用いる二つの関数をそれぞれ独立にフーリエ変換し, フーリエ変換結果をかけ合わせればよい, ということです.
もしくは, 畳み込みの計算はその逆方向フーリエ変換なので∫(f ∗ g)(x) =dτ f (τ )g(x − τ )= F −1 [f (k)g(k)](x)(2.12a)(2.12b)で計算できます.

2.2.3

畳み込みの例

ここでは畳み込みの例を考えてみます. ２つの関数

f (x) = x5 e−xg(x) = e−4x22の畳み込みを考えてみます.Eq. (2.8a) に代入してそのまま積分を計算すると,∫
h(x) =
f (x) · g(x − τ )dτ∫22=
τ 5 e−τ · e−4(x−τ ) dτ
√
4 2
1
π
x(1024x4 + 1600x2 + 375)e− 5 x=
3125 5(2.13)(2.14)(2.15a)(2.15b)(2.15c)となります.
一方で, それぞれの関数のフーリエ変換を掛け合わせて逆変換を行う Eq. (2.11), (2.12a)に従って求めてみたいと思います. Eq. (2.13) のそれぞれのフーリエ変換は,

iπ 3/22 2(4π 4 k 5 − 20π 2 k 3 + 15k)e−π k4
1 √ −π2 k2 /4
g(k) =
πe2

f (k) = −8

(2.16)(2.17)

2 連続の場合

2.2 畳み込み (Convolution)より 6 ,

h(k) = f (k)g(k)(2.18)

を逆変換して,∫ ∞dkh(k)ei2πkx(2.19a)

−∞
(
3/2

)
iπ
64ix(64x4 − 400x2 + 375) −4x2 /5 32ix(8x2 − 15) −4x2 /5 12ix −4x2 /5
√
√
= −
e
+
e
+√ e
(2.19b)43125 5π
25 5π
5π
√ (
)
2
1
π
=
x 1024x4 + 1600x2 + 375 e−4x /5(2.19c)3125 5と求められ, Eq. (2.15c) と一致することが確認できます.

6

f (x) の フ ー リ エ 変 換 https://www.wolframalpha.com/input?i=integral+from+-infty+to+

infty+of+x%5E5+exp%28-x%5E2%29*exp%28-i*2*pi*k*x%29g(x) のフーリエ変換 https://www.wolframalpha.com/input?i=integral+from+-infty+to+infty+of+
exp%28-4*x%5E2%29*exp%28-i*2*pi*k*x%29

93 離散の場合

3

離散の場合

3.1

離散フーリエ変換 (Discrete Fourier-Transform, DFT)
離散フーリエ変換とは, 位置・波数の両方が有限区間に制限された中で行われるフーリエ変換と考えて良いです.
つまり, 適当な離散的な基底関数 φn (x), (n = 0, 1, · · · ) を使用し, 離散的な点 x =xm , (m = 0, 1, · · · ) 上の値を用いて関数 f (xm ) を∑
f (xm ) =f (kn )φn (xm )

(3.1)

n

として書いたときに, f (xm ) と φn (xm ) が既知で, 係数 f (kn ) が未知の時, これを求めたいというのが離散フーリエ変換です. 本節の出発点は Eq. (2.1) とし, 順に離散的な表現に移行していきます. 順方向離散フーリエ変換は, ユニタリー性を保つために係数倍の自由度を指定する必要がありますので, 後程の項 3.1.3 で説明します.
本で離散フーリエ変換を知りたい方は, [2] などが良いかと思います.

3.1.1

位置を無限区間から有限区間へ制限 → 波数の刻み幅 ∆k の決定

まずは位置だけを区間 [xa , xb ], (xb > xa ) に制限して, その他の区間では [xa , xb ] が滑らかに繋がるように周期的に存在していると仮定します.
連続の場合の類推から, 離散フーリエ変換は基底関数

φn (x) = ei2πkn (x−xoﬀ )(3.2)

への射影となります 7 . ここで xoﬀ は基底関数を平行移動させる任意の値 (xoﬀ = 0 でも,

xoﬀ = xa でも OK), 波数 kn は離散的で,kn = n · ∆k, (n = 0, ±1, ±2, · · · )1∆k ≡xb − xa

(3.3)
(3.4)

だけが存在できます. ここで, ∆k は波数の間隔です. k が離散的であることと ∆k の値はどこから決まったかというと, 位置 [a, b] で周期的であるという条件から決められます. つまり, 区間の端で滑らかにつながるための周期境界条件


 φn (xa ) = φn (xb )dφn (x)dφn (x)= dx
dx x=xb
x=xa
7

xoﬀ の oﬀ は oﬀset の意味で使用しています.

10(3.5a)(3.5b)

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)を満たすような波数 k = kn しか存在しないと考える訳です. この答えが Eq. (3.3) になる訳です.
このままでも良いのですが, 位置だけ制限してなぜ波数は無限のままにするのか？という疑問が出てきます 8 . なので波数も有限の区間で考えましょう.

3.1.2

波数を無限区間から有限区間へ制限 → 位置の刻み幅 ∆x の決定基底関数 (3.2) に波数の最小値・最大値を制限した場合, どのような制限が掛かるか考
えてみましょう. Eq. (3.3) の制限から, 波数 k の最小値を kmin , 最大値を kmax とします.
式で表せば,

kmin = Nmin ∆k(3.6a)

kmax = Nmax ∆k(3.6b)

です. ここで, Nmin , Nmax はそれぞれの波数に対応するインデックスで Nmax > Nmin です.
波数 k の取り得る範囲とその範囲の大きさ krange は

k = [kmin , kmax ] = [Nmin ∆k, Nmax ∆k](3.7)

krange ≡ kmax − kmin = (Nmax − Nmin )∆k(3.8)となります.
波数の区間を制限すると, その範囲外の波数を持つ関数は範囲内に強制的に収められてしまいます. つまり, 範囲外の波数を持つ基底関数と範囲内の波数を持つ基底関数とが区別されなくなることを意味します. 式で表せば,

ei2πkn (x−xoﬀ ) = ei2π(kn +K)(x−xoﬀ ) ,(3.9)

ただし, x = [xa , xb ] で, kn + K は Eq. (3.7) の範囲外となってしまう定数 K が存在する, ということです. この K の値はいつも Eq. (3.7) の範囲外であるという条件から,

K = (Nmax − Nmin )∆k8(3.10)

対称性が悪い, とも言いますね. ある関数は位置で考え無ければならない, という制限はどこにもないのです. 多くの人が扱いやすいと感じるだけだから位置に注目しているだけであって, ある関数を波数で考えてもいいじゃないかという人もいるのです. どちらも優劣が無いので, 片方を制限したらもう片方も制限しない理由が特に無いのです. 逆に, わざと片方を連続のまま扱うことも可能で, 実際に量子力学の分野では利用されてたりします.

11位置 における値 ( )3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

( )

( )
( )

(

−1 )

3

2

( )
1

(

( )

)

( )
0

0

Δ

0

1

2

3

−1

−2

Δ
Δ

Δ

Δ

Δ

に
掛かる係数 ( )

( )
(
(
(

m  +2

m 

m 

+1

)

)

)

Δ
(

m

−1 )

(

m

)

(

m 

)

m 

m 

m

+1

m

−1

−1

m  +2Δ

Figure 1: 項 3.1.1, 3.1.2 で定義した分点の位置, 関係式等の図示.

となることは想像に難くないでしょう. 区間より大きくとって K > (Nmax − Nmin )∆k としても良いですが, 大きくとった際のメリットは特にないので通常は等号が成り立つときを K とします.
さて, Eq. (3.10) のよう K を決定すると, 具体的に kn のインデックス n の最大値を決めることができます. 安直に考えれば, n = Nmin , Nmin + 1, · · · , Nmax − 1, Nmax ですが, 端の Nmin と Nmax に置ける関数値は周期境界条件 (3.9) により

f (Nmin ∆k) = f (Nmax ∆k + K)

(3.11)

が成立しています. なので, 実際のインデックスは Nmin から Nmax − 1 までで打ち止めとなります. つまり, 波数はkn = n∆k + l · K,
(n = Nmin , Nmin + 1, · · · , Nmax − 1)

(3.12a)(3.12b)となるわけです. ここで, 波数は K の整数倍 l, (l = 0, ±1, · · · ) だけずれても良いので, それを考慮しています.
また Eq. (3.9) より

ei2πK(x−xoﬀ ) = 1 → K(x − xoﬀ ) = m, (m は整数)

12

(3.13)

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)の関係式がいつでも成立していなければなりません.

K, m は既に決まっているため,

Eq. (3.13) は x に対して制限を加える条件になっています. つまり, x − xoﬀ はいつでも 1/K(≡ ∆x) の整数倍になっていなければならない, という条件です. 式で表せば,

xm = m∆x + xoﬀ , (m = 0, 1, · · · , Nmax − Nmin − 1)
xb − xa1=∆x =K
Nmax − Nmin

(3.14a)(3.14b)となります. ここで, Nmax −Nmin −1 の ”−1” が入っているのは, m = 0 と m = Nmax −Nminが周期境界条件からいつも同じ値になるため, 除いても問題ないためです. また K を

Eq. (3.10) のように決定したことにより, xoﬀ = xa の制限が付きました. m = 0 の時に
x0 = xoﬀ となり, x0 は範囲の左端である xa に等しいからです.
以上で求めた関係について, Fig. 1 に示しましたので参考にしてください.

13

3 離散の場合

3.1.3

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

順方向/逆方向離散フーリエ変換の表記

フーリエ級数 f (kn ) による f (xm ) の表現は Eq. (3.1) のように書けます. そこでは具体的な点数などは決めていませんでしたが, ここまで来てようやくそれらを指定することができました.
逆方向離散フーリエ変換は Eq. (3.1) の通りで, 分点の位置や数は Eq. (3.12b) に従います. つまり,f (xm ) =

Nmax
∑−1

f (kn )ei2πkn (xm −xa )(3.15a)n=Nminとなります. 先に述べた xoﬀ は xa に等しいという事実を用いました. 順方向離散フーリエ変換を求めるためには, 連続フーリエ変換を求めた際に使用した Eq. (2.5) を真似すればよいです. Eq. (3.15) の両辺から e−i2πkn′ (xm −xa ) を掛けて m に対して和を取ってみます. す
ると,

N
−1
∑

e−i2πkn′ (xm −xa ) · f (xm ) = N · f (kn′ )(3.16)

m=0となり, 順方向離散フーリエ変換が導かれます. 計算過程は Appendix A.2 を参照してください.
以上から, 順方向/逆方向離散フーリエ変換は

N −1

1 ∑



f (kn ) =f (xm )e−i2πkn xm


N m=0
Nmax

∑−1



f (kn )ei2πkn xm
 f (xm ) =(3.17a)(3.17b)n=Nminとなります.
連続フーリエ変換 (2.1) とは異なり係数 1/N が存在し, 非常に, 非常に気持ち悪いです.
ですが, この N1 が存在しないとフーリエ変換のユニタリー性が保たれなくなります. 式で書くなら次の通りです.
−1

F ′ F ′ f (x) = N f (x)(3.18)

ここで, F ′ , F ′ −1 は Eq. (3.17) で 1/N が存在しない順方向, 逆方向離散フーリエ変換です.
これの右辺に入る N を出てこないようにするために両辺に 1/N を掛け, この定数が離散
フーリエ変換の定義に入り込んできます. そのため, 1/N は必ず順方向に入る訳ではなく,逆方向の方に入っても良いですし, 両方に入っても良いです.
この係数 1/N と指数関数の肩に 2π を含めるか否かの派閥のせいで, 本によって離散フーリエ変換の定義がバラバラです. 同じ著者でも気を付けてください. いくつか例を挙げてみましょう.

14

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

例 a) 対称性を重視する形式
例えば数学的に最もきれいな形にしようとすれば対称性を考えて

N −1

1 ∑


√

f (xm )e−i2πkn xm
f
(k
)
=
n


N m=0
N ∑−1


1 max


√
f (kn )ei2πkn xm
 f (xm ) = N
n=Nmin(3.19a)(3.19b)が良いでしょう. ですがどうもなかなか好まれないようです. 理由は様々ですが, ルー
トを書くのが面倒, プログラムの出力を見た際に分かりにくいというのがあります.
フーリエ変換のテストとして, f (xm ) = ei2πk0 x を実行してみると, Eq. (3.19a) の解√√は f (k0 ) = N となります. 例えば N = 5000 とかで実行したら 5000 の答えは計√算しなきゃいけません. ですが 1/ N が無ければ, そのまま f (k0 ) = 5000 が値とし√て出て嬉しいため, 1/ N はなかなか採用され無いのかな…と思います.
例 b) 連続フーリエ変換との対応を重視する形式
連続フーリエ変換の形と見比べると, 以下の形式でも良いことが分かります.

N
−1
∑




f (kn ) =f (xm )e−i2πkn xm ∆x(3.20a)







 f (xm ) =


m=0
Nmax
∑−1

f (kn )ei2πkn xm ∆k(3.20b)n=Nminこれは, フーリエ変換の積分を区分求積のように考えて離散フーリエ変換にしているように見えます. この書き方は連続フーリエ変換の答えがそのまま離散フーリエ変換の答えとなり, 係数倍の違いも生じません. つまり, ∆x∆k = 1/N が成立する
ため, これも意味がある分け方です. そのため連続フーリエ変換の純粋な離散表現を行いたい場合はこれが良いでしょう. しかし, マイナーです. 好んで使う人はなか
なかいない気がします (著者はこの表記が大好きです).
例 c) 可読性を重視する形式
もっともよく使用されるのが逆方向に N1 を入れる方法です.

N
−1
∑




f (kn ) =f (xm )e−i2πkn xm


m=0
Nmax
∑−1



1



 f (xm ) = N

f (kn )ei2πkn xm(3.21a)(3.21b)n=Nmin順方向をできるだけ簡単に表記しつつ, 逆方向に全てを押し付けた形ですね. 次の項以降でもこれを考えて行こうと思います.

15

3 離散の場合

3.1.4

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

離散フーリエ変換の簡略化

離散フーリエ変換に制限を加えることで, 表記がもう少し簡単になります. まずはこれまでの結果をまとめてみましょう. 前節 3.1.1, 3.1.2 の結果より, 離散フーリエ変換に関係する式は次の Eq. (3.22) から Eq. (3.27) の通りとなります.

1. 離散フーリエ変換


N
−1
∑




f (kn ) =f (xm )e−i2πkn xm


m=0
Nmax
∑−1



1


f
(x
m) =


N

f (kn )ei2πkn xm(3.22a)(3.22b)n=Nmin2. 位置 区間 [xa , xb ]xm = m∆x + xa , (m = 0, 1, · · · , N − 1)1
xb − xa∆x ==KN
N ≡ Nmax − Nmin

(3.23a)(3.23b)(3.23c)

3. 波数 区間 [kmin , kmax ], K = kmax − kmin = N ∆kkn = n∆k + lK, (n = Nmin , Nmin + 1, · · · , Nmax − 1)
∆k =

K
1
=
N
xb − xa(3.24a)(3.24b)4. 関数の周期境界条件
f (xm ± X) = f (xm ),(for all m)(3.25a)f (kn ± K) = f (kn ),(for all n)(3.25b)5. 関係式∆x∆k =

1, XK = NN(3.26)

6. 基底関数φn (x) = ei2πkn (x−xa )16(3.27)

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

離散フーリエ変換の指数関数部を詳細に見ていきましょう. つまり,

ei2πkn xm(3.28)に注目します. 具体的な表記 Eq. (3.23), (3.24) を代入していくと,

ei2πkn xm = ei2π(n∆k+lK)(m∆x+xa )(3.29a)= ei2πnm∆k∆x · ei2πlmK∆x · ei2π(n∆k+lK)xa(3.29b)= ei2πnm/N · ei2π(n∆k+lK)xa(3.29c)となります. ここで, ∆x∆k = 1/N, ∆x = 1/K を用いました. もし xa = 0 であれば非常に簡単になりますので, xa を 0 にとりましょう. つまり, 基底関数はいつも

φn (x) = ei2πkn xであり, 離散フーリエ変換は

N
−1
∑




f (kn ) =f (xm )e−i2πnm/N(3.30)(3.31a)

m=0

N
+N −1


1 min∑
f(x)=
f (kn )ei2πnm/N
m


N n=N

(3.31b)minと書くことができます. しかしながら, 位置の分点は

xm = m∆x, (m = 0, 1, · · · , N − 1)

(3.32)に変更されました. つまり, x 軸の範囲はもはや [xa , xb ] と任意の範囲ではなくなり, [0, N ∆x]の範囲でいつも記述されていなければならないことを意味しています. 簡単にいうなれば,

xm=0 = 0 を満たす決め方しか許されない代わりに簡単な表現 (3.31) を手に入れた, ということです. 良いところは, 分点の位置に左右されないことですね.∆x や ∆k など x, k の情報に左右されないのです.
また, 波数についても範囲を決めましょう. 波数は範囲が K = kmax − kmin であること
以外制限がなく, 決めた範囲以外は K で周期的なので, これまでの議論から一意に決めることはできません.

Nmin の決め方や, xm=0 ̸= 0 とする取り方にしたい場合にはフーリエ変換の表記が変わりますので, それらについて一通り紹介します.

17

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

• Nmin = 0 にとる場合
まずは最も多用される, Nmin = 0 とする場合を考えます. 何が便利かといえば, 逆離散フーリエ変換の式で添え字 n が 0 から N − 1 として表現でき, すっきりした表記になることが便利です. 波数は Nmin = 0 に決め, 以下の通りになります.

k = [lK, (l + 1)K],(3.33a)kn = n∆k + lK, n = 0, · · · , N − 1,

(3.33b)

K = N ∆k, ∆k =

K
1
=
N
N ∆x

(3.33c)となります. 多くの場合で l = 0 が選ばれます. 波数の範囲は Nmin によって自動的に決まるので, 離散フーリエ変換は,

N
−1

∑



f (kn ) =f (xm )e−i2πnm/N


m=0
N
−1
∑



1


f (kn )ei2πnm/N f (xm ) = N
n=0(3.34a)(3.34b)と書くことができます.x, k の範囲はそれぞれ x = [0, X], k = [0, K] です. 範囲は

xm=0 = 0, kn=0 = 0 である必要があります. 波数 k については特に考えずに周期境界条件を適用しても構いません 9 . つまり, k の範囲を k = [−K/2, K/2] で考えたい場合には添え字は kn=0 = 0 を満たすように,

k = [−K/2, K/2], n∆k,
n = 0, · · · , N/2 − 1,

N = 偶数の時 n∆k − K, n = N/2, · · · , N − 1,


kn =

 n∆k,
n = 0, · · · , (N − 1)/2,



N
=
奇数の時


 n∆k − K, n = (N − 1)/2 + 1, · · · , N − 1,

(3.35a)(3.35b)とすればよいです. 離散フーリエ変換は Eq. (3.34) を使用できます.

9

位置については, 不連続関数を取った場合に不連続点の位置が変わってしまいます. そのため, 考えて周期境界条件を適用しなければなりません. 詳細は, 次の項目位置の始まりを変更したい場合を参照してくだ
さい.

18

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

• 位置の左右対称に変更したい場合
最後に位置の始まりを左右対称にしたい場合を考えます. これは Eq. (3.32) で指摘したとおり, Eq. (3.29c) の後ろの指数関数がゼロにならなくなり (xa ̸= 0 となり),表現が複雑になってしまいます. 勝手に変更することはできません. それを承知のうえでも xm=0 = −X/2 のようにして, x = [−X/2, X/2] で考えたい場合を考えていきましょう.
これを解決するためには関数を −X/2 だけ平行移動してしまえばよいのです. つまり, 今考えたい新しい位置の分点 x′m は

x′ = [−X/2, X/2],(3.36a)x′m = m∆x − X/2, m = 0, · · · , N − 1,

(3.36b)としたいです.Eq. (3.22) に立ち戻ると, 今の状況は xa = −X/2 に相当するので′

ei2πkn xm を考えていくと′

ei2πkn xm = ei2πnm/N · ei2π(n∆k+lK)(−X/2)(3.37a)= ei2πnm/N · ei2π(−nN ∆k∆x/2−lXK/2)(3.37b)= ei2πnm/N · e−iπ(n+lN )(3.37c)= ei2πnm/N · (−1)n+lN(3.37d)となります. 以上から, 区間 [−X/2, X/2] で定義された関数の離散フーリエ変換はEq. (3.34) に代入すれば, 離散フーリエ変換は

N
−1

∑

 (−1)n+lN f (k ) =f (x′m )e−i2πnm/N
n









m=0
N
−1
∑

(3.38a)1f (x′m ) =(−1)n+lN f (kn )ei2πnm/N
N n=0となります. ここで, x は Eq. (3.36) に従い決定されます. 波数の範囲は Eq. (3.33) を考えて k = [0, K] としても良いですし, Eq. (3.35) に従って k = [−K/2, K/2] としても構いません (l = 0 を取った場合). 波数の周期境界条件を利用して範囲をずらす場
合は, 位置の場合と違って離散フーリエ変換の係数に影響を及ぼさないからです. 別の言い方では, 位置の起点 xa は基底関数の起点と関係があったので簡単にはいきま
せんが, 波数の起点は基底関数に関係ないです.
以上で, 位置が [0, X] もしくは [−X/2, X/2] で使えることを見てきましたが, もっと
一般的な範囲では見ていないので, それを考えてみましょう. これを考える意義は例えば,

19

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)x = 0 の近傍ではなくて遠い地点に波形の大部分がある場合, 離散フーリエ変換 (3.31) は使用できないのではないか？という問題に応えるためです. 既にみたように, Eq. (3.22) のxa を変更すれば良いのです.
これは平行移動の処理をしてしまえばよく, 平行移動した分だけを別で与えれば良いのだ, ということで解決することができます.
このアイデアを具体的に数式化しましょう. 本当にフーリエ変換したい波形 f (x) が定義域 [xa , xa + X] で定義されているとします. つまり,

x′ = [xa , xa + X],(3.39a)x′m = m∆x + xa , m = 0, · · · , N − 1,

(3.39b)′というわけです. ei2πkn xm は′

ei2πkn xm = ei2πnm/N · ei2πkn xa(3.40)

以上に進めることができません. よって, 離散フーリエ変換は

N
−1

∑
 i2πkn xaef (kn ) =f (x′m )e−i2πnm/N










f (x′m ) =

m=0
N
−1[
∑

]
1
i2πkn xaef (kn ) ei2πnm/N
N n=0(3.41a)の形となります. 通常の数値計算などでは, Eq. (3.41a) の右辺の形で実装されていることが多いです. よって, 本来のフーリエ変換結果である f (k) が欲しい場合は, 左辺に現れる余計な位相をキャンセルするようにこの因子を掛けなければなりません.

20

3 離散の場合

3.1.5

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

本稿の離散フーリエ変換の定義

本稿の離散フーリエ変換は Eq. (3.34) に従って以下の通りに定義します.

• 離散フーリエ変換


N
−1

∑


 f (kn ) =f (xm )e−i2πnm/N


m=0
N
−1
∑



1



 f (xm ) = N(3.42a)f (kn )ei2πnm/N(3.42b)xm = m∆x, (m = 0, 1, · · · , N − 1)

(3.43a)n=0• 位置 区間 [0, X]

∆x =

X
1
=
N
K

(3.43b)• 波数 区間 [lK, (l + 1)K],kn = n∆k + lK, (n = 0, 1, · · · , N − 1)
K
1
= ,
N
X
N
K = N ∆k =X
∆k =(3.44a)(3.44b)(3.44c)l = 0, ±1, ±2, · · · は任意の整数で, 典型的には l = 0 がとられ, 本稿でも注記が無い限り l = 0 とします. この l の値によってオーバーサンプリングやダウンサンプリングとして解釈されます.

• 関数の周期境界条件
f (xm ± X) = f (xm ),(for all m)(3.45a)f (kn ± K) = f (kn ),(for all n)(3.45b)• 関係式∆x∆k =

1, XK = NN(3.46)• 基底関数φn (x) = ei2πkn x ,(3.47a)φn (xm ) = ei2πnm/N ,(3.47b)特に x = xm の分点上では,となります.

21

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

任意の区間におけるフーリエ変換を知りたい場合は項 3.39 で示した通りの方法と注意で実装できます. また, 今後多くの場合で採用されている l = 0, N = (偶数) に絞り考えて
いきます.

22

3 離散の場合

3.1.6

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

離散フーリエ変換の例

項 2.1.3 と同じように, 同じ関数を用いて離散フーリエ変換を実施してみましょう.

f (x) = cos(2π3x)(3.48)

2π= 1/3 を持つ波) を離散フーリエ変換してみまという, 波数 3 を持つ波 (周期 Xp = 2π×3

す. ここでは, 離散フーリエ変換の区間 Eq. (3.43) と f (x) の周期が一致する場合と, 一致
しない場合の２つの例を考えてみます.

• 離散フーリエ変換の区間 Eq. (3.43) と f (x) の周期が一致する場合
離散フーリエ変換の区間を x = [0, qXp ] と波形の周期の整数倍として考えてみます.
ここで, q は正の整数で q = 1, 2, · · · を取るものとします. 離散フーリエ変換における特徴的な量は以下のように決まります.

qXp
1, ∆k =
,
N
qXp
[
]
N
N
K=
→ k = 0,
qXp
qXp
∆x =

(3.49)(3.50)

まずは順方向離散フーリエ変換をしてみます.Eq. (3.42a) に代入して,

f (kn ) =
=

=

=
=

N
−1
∑
m=0
N
−1
∑
m=0
N
−1
∑

f (xm )e−i2πnm/N
cos(2π3 · mcos(2π

m=0
N
−1 [
∑

qXp −i2πnm/N)eN

qm −i2πnm/N)e
N(3.51a)(3.51b)(3.51c)

]
(q+n)
(q−n)
1
ei2πm N + e−i2πm N
2 m=0(3.51d)N(δn,q + δn,−q )2(3.51e)となります. ここで, δi,j はクロネッカーのデルタで 1, (i = j)
δi,j =
 0, (i ̸= j)(3.52)

です.q > 0 であることから, Eq. (3.55) の右辺の δn,−q で n = −q になることはないと思われますが, 波数に対する周期境界条件 f (kn ) = f (kn ± K) があることから,

23

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

δn,−q+N とするべきで, 値を持つことになります. つまり, 現在の定義で順方向離散フーリエ変換の結果は

f (kn ) =N
(δn,q + δn,−q+N ) , (n = 0, 1, · · · , N − 1)
2

(3.53)となります. ここで, 等比級数の総和を用いて
N
−1
∑

′

ei2πm(n−n )/N = N δn,n′(3.54)

m=0を用いました.

Eq. (3.53) より, n = q, N − q の位置にピークが現れることが分かります. これはつまり, 波数としては

kn = n∆k =1 q∆k = q qX = 3
p

1


=K −3
 (N − q)∆k = N ∆k − q
qXp(3.55)となり, 波数 3 と K − 3 にピークが出ることが分かります. 周期境界条件により, 波
数に K の整数倍を足したり引いたりしても良いので, 分かりやすい表現では, 波数の成分 ±3 を持つことが分かります. この 3 は元々の波形の波数が 3 だから, ということで期待する通りの結果であることが分かります.
最後に連続フーリエ変換との比較をしましょう. 連続フーリエ変換の結果は Eq. (2.7)の通り, デルタ関数を用いて

1
f (k) = [δ(k − 3) + δ(k + 3)]
2(3.56)となります. それに対し, 離散フーリエ変換では Eq. (3.53) の通りで

f (kn ) =N(δkn ,3 + δkn ,−3 )2(3.57)

です. 両者を比較するとデルタ関数がクロネッカーのデルタに対応していることが分かります. 実際には比較するためには, 離散フーリエ変換で N1 に押し込めた ∆x∆kのうち, ∆x を掛けなければなりません (Eq. (3.20a) を参照). つまり,

δ(k − 3) → N δkn ,3 ∆xに対応することが分かります.

24(3.58)

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)

• 離散フーリエ変換の区間 Eq. (3.43) と f (x) の周期が一致しない場合
続いて, 離散フーリエ変換の区間を x = [0, X] として, 任意の範囲で考えてみます.
特に, X が周期の無理数倍と仮定して, 必ず周期境界条件を満たさない場合を仮定します 10 .
離散フーリエ変換における特徴的な量は以下のように決まります.

∆x =
K=

X
1
, ∆k = ,N
X[

NN→ [0, K] = 0,X
X]

(3.59)(3.60)

順方向離散フーリエ変換をしてみます.Eq. (3.42a) に代入して,

f (kn ) =
=

N
−1
∑
m=0
N
−1
∑

f (xm )e−i2πnm/N
cos(2π3 · m

m=0
N
−1 [
∑

X −i2πnm/N)eN](3X−n)(3X+n)
1i2πm N−i2πm N
=
e
+e
2 m=0
( i2π(3X−n))e−11e−i2π(3X+n) − 1=+
2 ei2π(3X−n)/N − 1 e−i2π(3X+n)/N − 1()
1ei2π3X − 1
e−i2π3X − 1=+
2 e−i2π(kn −3)∆x − 1 e−i2π(kn +3−K)∆x − 1(3.61a)(3.61b)(3.61c)(3.61d)(3.61e)となります. 注記しますが, X は f (x) の周期の無理数倍であることから, 離散フーリ
エ変換で参照される波数 kn = n∆k = Xn が 3 の整数倍に等しくなることはないため,

Eq. (3.61e) の分母が発散することはなく, Eq. (3.61c) から Eq. (3.61d) の計算が可能になっています. また, Eq. (3.61e) の右辺第二項の分母に突如現れた K は有ってもなくても影響しませんが, 分母の発散が分かりやすいように, また定義域に収まるよ
うに追加しました.
では結果 (3.61e) が何を示すのか考えていきましょう. 右辺の分母に着目すると, kn =3, K − 3 に近い場合にゼロに近づくため, f (kn ) の値が大きくなることが予想されます. 以降, kn = 3, K − 3 というのは面倒なので, K だけずらして kn = ±3 と呼びます.
以上から, f (kn ) は kn = ±3 周りにある程度の幅を持ったピークを持つ関数となる
ことが分かります. 幅を持つということは, k = ±3 以外の成分が出ていることを意
10

周期境界条件を満たさない条件なら半整数倍, もしくは有理数倍にすればよいですが, そうしない理由について Appendix B をご覧ください.

25

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)DFT result of f(x)=cos(2π3x), Analytical solution
with x=[0,2exp(1)], N=200, K≈36.8
100|f(k)|
Re(f(k))Im(f(k))

f(k)

50
0
-50
-100
0

3

6

9

12

15

18
k

21

24

27

30

33

36

Figure 2: Eq. (3.48) の DFT 結果 (2) の図示. X = [0, 2 exp(1)], N = 200, K = 1/∆x ≈ 36.8 として計算味しています. 元々の波形 (3.48) には k = ±3 以外の成分が含まれていなかったのに, どうしてそれ以外の波数が出てくるのでしょうか. これは周期境界条件を満たした波形を離散フーリエ変換していないことが原因です. 離散フーリエ変換の基底関数が周期境界条件を満たしているため, そこから構成される f (x) も周期境界条件(3.5) を満たしていなければなりません. そうでなければ, 離散フーリエ変換は適用するべきではなく, Eq. (3.61e) は適用するべきではないときの結果となります.
しかし, この適用するべきではないという事実は理不尽です. なぜなら, 多くの場合
ではその関数の波数を知りたいために離散フーリエ変換を行うのですが, その離散フーリエ変換を行うためにはその関数の波数を知っていなければならないというのは理不尽極まりないのです. もし関数の波数を知らない場合に, 正確な離散フーリ
エ変換を行いたければ, 離散フーリエ変換の範囲 X は有理数で表現できなくなる値
まで大きくとらなければなりません. つまり, X → ∞ として無限大の区間で行わなければならず, これはもはや離散フーリエ変換ではなくなります.
多くの場合では有限に収めなければならないため, 妥協します. つまり, 正確ではな
いけど Eq. (3.61e) はほとんど正解に近いはずだと考える訳です. 図で Eq. (3.61e) を書けば, Fig. 2 のようになります. k = 3 と, k ≈ 33.8 に幅を持つピークが現れていることが分かります. もしも k の範囲を正負で考えて K だけずらせば, k = ±3 となります. これは, もともと考えたかった波形が持つピークの位置となります. ですが,幅を持って広がるため, Eq. (3.55) のようにその周波数だけピークになる, ということにはなりません.

26

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)
function for different range f(x)=cos(2π3x)with X=2exp(1), x=[0,X] and x=[-X/2,X/2]

1.5

x=[0,X]

1

f(x)

0.5
0
-0.5
-1
-1.5
-6
1.5

-5

-4

-3

-2

-1

0
x

1

2

3

4

5

6

-4

-3

-2

-1

0
x

1

2

3

4

5

6

x=[-X/2,X/2]

1

f(x)

0.5
0
-0.5
-1
-1.5
-6

-5

Figure 3: 周期関数の異なる範囲の離散フーリエ変換を考えた場合. (上)x = [0, X], (下)x = [−X/2, X/2].
色ごとに X の範囲を示す.

• 周期が一致しないかつ左右対称の区間を選んだ場合

今までは位置の範囲が x = [0, X] でしたが, x = [−X/2, X/2] の場合を考えてみましょう. この２つは同じ結論にはなりません. その理由は, 不連続点が x = 0 に来る
か,x = X/2 に来るかの違いがあるためです.Fig. 3 にその比較を示しました. この違いによって, 離散フーリエ変換が認識する関数の周期に差が生じます. この場合のフーリエ変換は, Eq. (3.38a) でみた形で計算しなければなりません. つまり, 離散

27

3 離散の場合

3.1 離散フーリエ変換 (Discrete Fourier-Transform, DFT)DFT result of f(x)=cos(2π3x), Analytical solution
with x=[-exp(1),exp(1)], N=200, K≈36.8
100|f(k)|
Re(f(k))Im(f(k))

f(k)

50
0
-50
-100
0

Figure 4:

3

6

9

12

15

18
k

21

24

27

30

33

Eq. (3.48) を左右対称の区間で離散フーリエ変換した解析解 (3.63) の図示.

36

X =

[− exp(1), exp(1)], N = 200, K = 1/∆x ≈ 36.8 として計算

フーリエ変換の結果は l = 0 を考えて,(−1)n f (kn ) =
=

N
−1
∑
m=0
N
−1
∑

f (xm )e−i2πnm/N(3.62a)cos(2π3 · (m∆x − X/2))e−i2πnm/N(3.62b)

m=0

1
=
2

(
)
ei2π3X − 1
e−i2π3X − 1
−iπ3X
iπ3X
e
+e
(3.62c)e−i2π(kn −3)∆x − 1e−i2π(kn +3−K)∆x − 1

より,

(
)
1
ei2π3X − 1
e−i2π3X − 1
−iπ3X
iπ3X
f (kn ) = (−1) · · e+e(3.63)
2e−i2π(kn −3)∆x − 1e−i2π(kn +3−K)∆x − 1nとなります.Eq. (3.61e) と比べて (−1)n と, 指数関数が追加された形になります. 実
際に Eq. (3.63) を図示すると, Fig. 4 のようになります.

28

3 離散の場合

3.2

3.2 離散畳み込み (Convolution)

離散畳み込み (Convolution)
節 2.2 で考えた畳み込みを離散的に考えてみましょう. 離散畳み込みとは, 関数 f (x), g(x)の値が, 同じ区間で離散的かつ等間隔の離散点にだけ値が与えられているときに以下の和を計算することです.

∑(f ∗ g)d (xn ) ≡f (xl )g(xn − xl )(3.64a)lf (x), g(x) の定義域は同じ範囲です. 一点注意しておかなければならない点は, 連続の場
合に定義した畳み込み (2.8a) をそのまま離散的にしたのではないことです.Eq. (2.8a) をそのまま離散化するならば, 例えば(f ∗ g)(xn ) =≈∫ ∞f (τ )g(xn − τ )dτ

(3.65a)f (xp )g(xn − xp )∆x(3.65b)−∞∑p= (f ∗ g)d (xn ) · ∆x

(3.65c)となります. ですが, 多くの場合では Eq. (3.64a) のように ∆x を消去したものを離散畳み込みと呼んでいるので, 本稿もそれに倣います.
さて, Eq. (3.64a) で添え字の範囲などを明確にを指定しなかったのには理由があります. これは, Eq. (3.64a) の定義において xn − xl が定義域から飛び出してしまう場合が存在するためです. この飛び出した範囲の g(x) をどのように処理するかで, 畳み込みの結果が変わります. 主な決め方は２つあり, 関数を周期的と考える『循環畳み込み』, もしくは定義域外はゼロと考える『直線畳み込み』です.

• 循環畳み込み
循環畳み込みは定義外の領域を周期的であるとする方法です. つまり,

g(xn ± X) = g(xn ), (for all n)(3.66)

として考えます (f についても同様です). ここで, X は定義域の大きさです.

• 直線畳み込み
直線畳み込みは定義外の領域をゼロであるとする方法です. つまり,
 g(xn ), (xn は定義域内)
g(xn ) = 0,(otherwise)として考えます.

29(3.67)

3 離散の場合

3.2 離散畳み込み (Convolution)

本稿では特に循環畳み込みを考えていきます. これは離散フーリエ変換を利用して畳み込みを計算したいという意図があるため, 同じ周期境界条件を持つものを仮定したいた
めです.

3.2.1

離散畳み込みの定義

本稿における離散畳み込みは循環畳み込みを想定し, 以下のように定義します.

(f ∗ g)d (xm ) ≡N−1∑

f (xp )g(xm − xp )(3.68a)p=0

f (xm ± X) = f (xm )(3.68b)g(xm ± X) = g(xm )(3.68c)

ここで, X は定義域の区間であり, 位置の分点と共に以下の通り定義します.

x = [0, X](3.69a)xm = m∆x, (m = 0, 1, 2 · · · , N − 1)

(3.69b)∆x =X
N(3.69c)もし, 別の定義域で定義された波形を考えたい場合は, 平行移動して [0, X] の範囲に移動させて計算し, その後全て計算し終えたら平行移動した分だけ戻すことで実行できます.

3.2.2

離散畳み込みと離散フーリエ変換

項 2.2.2 の連続の場合で考えたように, 離散畳み込みを離散フーリエ変換された関数で考えてみましょう. 離散フーリエ変換の定義は Eq. (3.42a), (3.42b) に従います. 循環畳み込みを考えたいので, 関数は区間で周期的であるから関数 f, g は

1 ∑f (xm ) =f (kn )ei2πnm/N
N n=0(3.70a)1 ∑g(xm ) =g(kn )ei2πnm/N
N n=0(3.70b)

N −1

N −1

30

3 離散の場合

3.2 離散畳み込み (Convolution)と書けています.Eq. (3.70a), (3.70b) を Eq. (3.68a) に代入すれば,(f ∗ g)d (xn ) =

N
−1
∑
p=0

=

N
−1
∑
p=0

f (xp )g(xm − xp )[1 ∑f (kn )ei2πnp/N
N n=0
N −1

(3.71a)][1 ∑′g(kn′ )ei2πn (m−p)/N
N n′ =0
N −1](3.71b)

N
−1
N −1
N
−1
∑
∑
1 ∑
′
i2πn′ m/N=
ei2π(n−n )p/N
f (kn )
g(kn′ )e
2
N n=0
p=0
n′ =0

(3.71c)

N −1
N
−1
∑
1 ∑
′
=
f (kn )
g(kn′ )ei2πn m/N N δn,n′
2
N n=0
n′ =0

(3.71d)1 ∑=f (kn )g(kn )ei2πnm/N
N n=0(3.71e)N −1と計算することができます. つまり, 関数 f (kn )g(kn ) の逆変換が離散畳み込みの結果を与える, といっています. 結局, 項 2.2.2 と同じ結論が得られるわけですね. 畳み込みの離散
フーリエ変換の結果なので, (f ∗ g)d (kn ) ≡ f (kn )g(kn ) と書きましょう. 実際に畳み込みを
フーリエ変換すると以下のようになります.

(f ∗ g)d (kn ) ≡ f (kn )g(kn ) = F[(f ∗ g)d (x)](kn )
[ −1
]
N
−1 N
∑
∑
=
f (xp )g(xm − xp ) e−i2πnm/N
m=0

=

N
−1
∑

=

−i2πnp/N

f (xp )e

=

N
−1
∑

g(xm − xp )e−i2πn(m−p)/N (3.72c)m=0−i2πnp/N

f (xp )e

N
−1
∑

′

g(xm′ )e−i2πnm /N(3.72d)

m′ =0

p=0
N
−1
∑

(3.72b)

p=0

p=0
N
−1
∑

(3.72a)f (xp )e−i2πnp/N · g(kn )(3.72e)p=0= f (kn )g(kn )(3.72f)

ここで, Eq. (3.72c) から Eq. (3.72d) に移る際に循環畳み込みの性質 (3.68c) を利用しています. これも項 2.2.2 と同様の結論が得られたことになります. 畳み込み積分 (3.68a)
自体を知りたければ, Eq. (3.72f) を逆変換することで得ることができます.

3.2.3

任意区間における離散畳み込みと離散フーリエ変換

任意区間における離散畳み込みと離散フーリエ変換の関係を考えてみます.

31

3 離散の場合

3.2 離散畳み込み (Convolution)x = 0 を中心とした関数なので x = [xa , xa + X] で定義されたフーリエ変換 Eq. (3.41a)を元に考えていきます.l = 0 として再掲すれば

N
−1

∑

i2πk
x

n
a

e
f (kn ) =f (x′m )e−i2πnm/N


m=0
N
−1[
∑

(3.73a)]1f (x′m ) =ei2πkn xa f (kn ) ei2πnm/N
N n=0








位置, 波数の位置はそれぞれ Eq. (3.74), (3.75) のように決めています.

x′ = [xa , xa + X],(3.74a)x′m = m∆x + xa , m = 0, · · · , N − 1,

(3.74b)k = [0, K],(3.75a)kn = n∆k, n = 0, · · · , N − 1,

(3.75b)xm − xp = (m − p)∆x を計算する際は xa = 0 と考えることに注意して(f ∗ g)d (x′m )

=

N
−1
∑
p=0

=

N
−1
∑
p=0

=

f (x′p )g(x′m − x′p )[(3.76a)N −1

1 ∑ i2πkn xaef (kn )ei2πnp/N
N n=0

][

N −1

′
1 ∑
g(kn′ )ei2πn (m−p)/NN ′](3.76b)

n =0

N −1
N
−1
N
−1
∑
∑
′
1 ∑ i2πkn xa
i2πn′ m/N
′ )e
e
f
(k
)
g(k
ei2π(n−n )p/N
n
n
2
N n=0
′
p=0

(3.76c)

N −1
N
−1
∑
′
1 ∑ i2πkn xaef(k)
g(kn′ )ei2πn m/N · N δn,n′
n
2
N n=0
′

(3.76d)N −1][]
1 ∑ −i2πkn xa [ i2πkn xaeef (kn ) ei2πkn xa g(kn ) ei2πnm/N
N n=0(3.76e)

n =0

=

n =0

=

となります.Eq. (3.76e) の中括弧[]

内は離散フーリエ変換 Eq. (3.73) の結果に対応します. つまり, フーリエ変換を通して畳み込みを計算する場合は, 位相因子 e−i2πkn xa を付与しなければならないことに注意してください. 連続畳み込みと対応付ける際は, さらに積分の刻み幅である ∆x を最後の結果に掛けることも忘れないようにしましょう.

3.2.4

離散畳み込みの例

ここでは離散畳み込みの例を考えてみます. 問題は連続の畳み込み (2.2.3) と同じように２つの関数f (x) = x5 e−xg(x) = e−4x

32

2

2

(3.77)(3.78)

3 離散の場合

3.2 離散畳み込み (Convolution)の離散畳み込み(f ∗ g)d (xm ) =∑f (xp )g(xm − xp )(3.79)

pを考えてみます.
今回考える関数は x = 0 を中心とした関数なので, x = [xa , xa + X] で定義されたフーリエ変換 Eq. (3.73) を元に考えていきます. ここで x′m と kn は, Eq. (3.74), (3.75) で決めます.
それぞれの連続フーリエ変換を離散化することで簡易的に求めましょう. この対応は,−xa , xa + X >> 1 のみで成り立ちます. 連続の結果は項 2.2.3 で得た結果をただ離散化す
ればよいので,iπ 3/212 2(4π 4 k 5 − 20π 2 k 3 + 15k)e−π k ·(3.80)4∆x
1 √ −π2 k2 /4 1
g(k) =πe·(3.81)
2∆x∫∑となります. ここで 1/∆x は連続から離散の表現に伴って生じる, dx →∆x への変形f (k) = −時に生じる ∆x を左辺に移した形です. Eq. (3.79) の解も同じように離散化すれば,(f ∗ g)d (xm ) =∑f (xp )g(xm − xp )

p

1
=
3125

√

)
π (
1
2
xm 1024x4m + 1600x2m + 375 e−4xm /5 ·
5
∆x

を得ます.

33

(3.82)(3.83)

4 プログラム

4

プログラム
本章では畳み込みを実行するプログラムを載せます. 言語は Fortran90 です.

Intel MKL を利用した離散 (高速) フーリエ変換

4.1
4.1.1

Intel MKL の離散フーリエ変換の定義

Intel MKL の離散フーリエ変換は, Eq. (3.42) で見た形で定義されています.
• 離散フーリエ変換


N
−1
∑




f
(k
)
=
f (xm )e−i2πnm/Nn(4.1a)

m=0

N
−1

∑



f
(x
)
=
f (kn )ei2πnm/Nm


(4.1b)n=0• 位置
x = [0, X](4.2a)

xm = m∆x, (m = 0, 1, · · · , N − 1)

(4.2b)

ここで, ∆x は以下で定義されます.

∆x =

X
1
=
N
K

(4.3)

• 波数
k = [0, K](4.4a)

kn = n∆k + lK, (n = 0, 1, · · · , N − 1)

(4.4b)

ここで, ∆k, K は以下で定義されます.

K
1
= ,
N
X
N
K = N ∆k =X
∆k =(4.5a)(4.5b)• 関数の周期境界条件
f (xm ± X) = f (xm ),(for all m)(4.6a)f (kn ± K) = f (kn ),(for all n)(4.6b)• 関係式∆x∆k =

1, XK = N
N

34

(4.7)

4 プログラム

4.1 Intel MKL を利用した離散 (高速) フーリエ変換• 基底関数φn (x) = ei2πkn x ,(4.8a)

φn (xm ) = ei2πnm/N ,(4.8b)

特に x = xm の分点上ではとなります.

35

4 プログラム

4.1.2

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

プログラム

Intel MKL を利用したプログラムを書いてみましょう. 例えば Program (1) のようになります. 分かりやすさを中心に書いていますので, 無駄があります. プログラム中の変数 Np が 2, 3, 5, 7 の乗数になると高速フーリエ変換として自動的に実行されます. MKL の詳しい使い方はマニュアル [1] をご覧ください.
Listing 1: Intel MKL を用いた Fortran90 のプログラム例
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

program main
implicit none
integer::mdouble precision::dx,Xrange
double precision,allocatable::x(:)
integer::ndouble precision::dk,Krange
double precision,allocatable::k(:)
integer::Npcomplex(kind(0d0)),allocatable::f(:)
complex(kind(0d0)),external::func
Np = 200 ! Number of DFT points
Xrange = 2d0*exp(1d0) ! x=[0,Xrange]! Initialize and allocation for DFT
allocate(x(0:Np-1),k(0:Np-1),f(0:Np-1))
x = 0d0
k = 0d0
f = dcmplx(0d0,0d0)! Define position x
dx = Xrange / Np
call dft_abscissa(Np, x, dx) ! x=[0, Xrange]! Calculate k abscissa
Krange = 1d0/dx
dk = Krange / Np
call dft_abscissa(Np, k, dk) ! k=[0, Krange]! Calculate function f at x=x(m)do m = 0,Np-1
f(m) = func(x(m))
enddo! Write down f(x)call dft_write("before_fx.d",Np,x,f)! Forward dft (f is overwrited)call dft(Np,f,"forward")! Write down f(k)call dft_write("fk.d",Np,k,f) ! output k=[0,Krange]
call dft_write_center0("fk_center0.d",Np,k,f) ! output k=[-Krange/2,Krange/2]! Backward dft (f is overwrited)call dft(Np,f,"backward")! Write down f(x)call dft_write("after_fx.d",Np,x,f)
stop
end program main
complex(kind(0d0)) function func(x)implicit none
double precision,intent(in)::x!
! Defined x=[0,X]!
double precision::pi=dacos(-1d0)
func=dcmplx(dcos(2d0*pi*3d0*x),0d0)!func=dcmplx(dsin(2d0*pi*5d0*x),0d0)!func=exp(dcmplx(0d0,2d0*pi*1*x))

36

4 プログラム

66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

returnend function func!---------------------------------subroutine dft_abscissa(N, w, dw)
implicit none
integer,intent(in)::N
double precision,intent(in)::dw
double precision,intent(out)::w(0:N-1)!
!w=[0, W]! w(0) = 0! w(N-1) = (N-1)*dw = W-dw!
integer::j
do j = 0,N-1
w(j) = j*dw
enddo
returnend subroutine dft_abscissa
include "/opt/mkl/include/mkl_dfti.f90"
subroutine dft(N,f,FB)
use MKL_DFTI
implicit none
integer,intent(in)::N
complex(kind(0d0)),intent(inout)::f(0:N-1)
character(*),intent(in)::FB!
!sikinote! author : sikino! date : 2022/06/11!
!n: Number of DFT points.
!f(i): value of data at i
!FB: "forward" -> Forward DFT! "backward" -> Backward DFT
integer::Status
TYPE(DFTI_DESCRIPTOR),POINTER::hand
Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,N)
Status = DftiSetValue(hand,DFTI_FORWARD_SCALE,1d0)
Status = DftiSetValue(hand,DFTI_BACKWARD_SCALE,1d0/dble(N))
Status = DftiCommitDescriptor(hand)if(trim(FB)=="forward")then
Status = DftiComputeForward(hand,f)elseif(trim(FB)=="backward")then
Status = DftiComputeBackward(hand,f)else
write(6,*)"DFT string different"
stop
endif
Status = DftiFreeDescriptor(hand)
returnend subroutine dft
subroutine dft_write(fname,N,w,f)
implicit none
character(*),intent(in)::fname
integer,intent(in)::N
double precision,intent(in)::w(0:N-1)
complex(kind(0d0)),intent(in)::f(0:N-1)! w = [0,W]! w(0) = 0! --> output range w=[0,W]
integer::jopen(21,file=trim(fname))
do j = 0,N-1
write(21,’(3e26.16e3)’)w(j),dble(f(j)),dimag(f(j))
enddoclose(21)return37

4 プログラム

147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

end subroutine dft_write
subroutine dft_write_center0(fname,N,w,f)
implicit none
character(*),intent(in)::fname
integer,intent(in)::N
double precision,intent(in)::w(0:N-1)
complex(kind(0d0)),intent(in)::f(0:N-1)! w = [0,W]! w(0) = 0! --> output range w=[-W/2,W/2]
integer::j,Nhalf
double precision::dw,Wrange
if(mod(N,2)==0)then
Nhalf = N/2
else
Nhalf = (N-1)/2
endif
dw = abs(w(1)-w(0))
Wrange = N*dw
open(21,file=trim(fname))
do j = Nhalf,N-1
write(21,’(3e26.16e3)’)w(j)-Wrange,dble(f(j)),dimag(f(j))
enddo
do j = 0,Nhalf-1
write(21,’(3e26.16e3)’)w(j),dble(f(j)),dimag(f(j))
enddoclose(21)returnend subroutine dft_write_center038

4 プログラム

4.1.3

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

MKL を利用した離散フーリエ変換の計算例

実際に離散フーリエ変換をプログラム上で実行してみます. 波の形は Eq. (2.1.3) やEq. (3.48) で解いた時と同じ,
f (x) = cos(2π3x)

(4.9)

とします.

• 離散フーリエ変換の区間と f (x) の周期が一致する場合
波の周期と x 軸の範囲 X が整数倍なので, この解析解は Eq. (3.55) より,

N(4.10)(δn,q + δn,−q )2)N(=δkn ,kq + δkn ,k−q(4.11)
2となります. 具体的に X = [0, 5], N = 200, K = 1/∆x = 40 として計算することを
f (kn ) =考えると q = 15 なので, ±q∆k = ±3 となります. よって,f (kn ) = 100 (δkn ,3 + δkn ,−3 )(4.12)となります. 数値的に DFT を実行した結果は Fig. 5 のようになります. 数値計算上
でも Eq. (4.12) の通りになっていることが確認できます.

• 離散フーリエ変換の区間と f (x) の周期が一致しない場合
Eq. (4.9) で DFT の区間が波の周期の無理数倍である時かつ, 区間が x = [0, X] の場
合, 解析解はf (kn ) =N−1∑f (xm )e−i2πnm/N

m=0

1
=
2

(

ei2π3X − 1
e−i2π3X − 1+
e−i2π(kn −3)∆x − 1 e−i2π(kn +3)∆x − 1(4.13a))(4.13b)となります. 数値的に DFT を実行した結果は Fig. 6 のようになります. 図中でフーリエ変換結果を比較している通り Eq. (4.13) と解析解が一致していることが分かります.
また, Eq. (4.9) で DFT の区間が波の周期の無理数倍である時かつ, 区間が x =[−X/2, X/2] の場合, 解析解は Eq. (3.63) で見た通り()
ei2π3X − 1
e−i2π3X − 1
−iπ3X
iπ3X
n 1
+e
f (kn ) = (−1) · · e(4.14)
2e−i2π(kn −3)∆x − 1e−i2π(kn +3−K)∆x − 1となります. 数値的に DFT を実行した結果は Fig. 7 のようになります. 図中でフーリエ変換結果を比較している通り Eq. (4.14) と解析解が一致していることが分かります.

39

4 プログラム

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

DFT result for f(x)=cos(2π3x)with X=5, x=[0,X], N=200, K=40.000
1(before DFT)

f(x)

0.5

Re(f(x))Im(f(x))

0
-0.5
-1
0

1

2

3

x

4

5

f(k)

100
50

(DFT result)0

Re(f(k))Im(f(k))

-50
-20

-15

-10

-5

0
k

5

10

15

20

1
(IDFT result)

f(x)

0.5

Re(f(x))Im(f(x))

0
-0.5
-1
0

1

2

3

x

4

5

Figure 5: Eq. (4.9) の DFT 結果. X = [0, 5], N = 200, K = 1/∆x = 40 として, プログラム (1) を用いて計算.

40

4 プログラム

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

DFT result for f(x)=cos(2π3x)with X=2exp(1), x=[0,X], N=200, K≈36.788
1(before DFT)

f(x)

0.5

Re(f(x))Im(f(x))

0
-0.5
-1
0

1

2

x

3

4

5

100
(DFT result)

f(k)

50
0

Re(f(k))Im(f(k))Re(fA(k))
Im(fA(k))

-50
-15

-10

-5

0
k

5

10

15

1
(IDFT result)

f(x)

0.5

Re(f(x))Im(f(x))

0
-0.5
-1
0

1

2

x

3

4

5

Figure 6: Eq. (4.9) の DFT 結果. X = [0, 2 exp(1)], N = 200, K = 1/∆x ≈ 36.8 として, プログラム (1)を用いて計算. 図中 fA (k) は, 解析解 (4.13) を表す.

41

4 プログラム

4.1 Intel MKL を利用した離散 (高速) フーリエ変換

DFT result for f(x)=cos(2π3x)with X=2exp(1), x=[-X/2,X/2], N=200, K≈36.788
1(before DFT)

f(x)

0.5

Re(f(x))Im(f(x))

0
-0.5
-1
-2

-1

0
x

1

2

100
(DFT result)

f(k)

50
0

Re(f(k))Im(f(k))Re(fA(k))
Im(fA(k))

-50
-15

-10

-5

0
k

5

10

15

1
(IDFT result)

f(x)

0.5

Re(f(x))Im(f(x))

0
-0.5
-1
-2

-1

0
x

1

2

Figure 7: Eq. (4.14) の DFT 結果. X = [− exp(1), exp(1)], N = 200, K = 1/∆x ≈ 36.8 として, プログラム (5) を用いて計算. 図中 fA (k) は, 解析解 (4.14) を表す.

42

4 プログラム

4.2 Intel MKL を利用した畳み込みIntel MKL を利用した畳み込み

4.2
4.2.1

Intel MKL を利用した畳み込みの計算例

ここでは項 3.2.4 で計算した具体的に畳み込みを計算します. つまり,

f (x) = x5 e−xg(x) = e−4x22(4.15)(4.16)の離散畳み込みhA (xm ) ≡ (f ∗ g)d (xm )∑=f (xp )g(xm − xp )
p

1
=
3125

√

)
π (
1
2
4
2
xm 1024xm + 1600xm + 375 e−4xm /5 ·5∆x

(4.17)(4.18)(4.19)

を考えます. 途中計算で現れる形も, 項 2.2.3 と同じように示せば

iπ 3/212 2(4π 4 k 5 − 20π 2 k 3 + 15k)e−π k ·4∆x1 √ −π2 k2 /4 1
gA (k) =πe·2∆xfA (k) = −となります.
離散フーリエ変換の結果は Fig. 8 となります.

43(4.20)(4.21)

4 プログラム

4.2 Intel MKL を利用した畳み込み

22Discrete convolution f(x)=x5e-x and g(x)=e-4xwith X=20, x=[-X/2,X/2], N=200, K≈10.000

f(x)

1
0.5

(before DFT)

0

f(x)
g(x)

-0.5
-1
-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10x
20
(DFT of f(x))
f(k)

10

Re(f(k))Im(f(k))Re(fA(k))
Im(fA(k))

0
-10
-20

-5

-4

-3

-2

-1

0
k

1

2

3

4

5

20
(DFT of g(x))
g(k)

10

Re(g(k))Im(g(k))Re(gA(k))
Im(gA(k))

0
-10
-20

-5

-4

-3

-2

-1

0
k

1

2

3

4

5

10
(IDFT of h(x))h(x)

5Re(h(x))Im(h(x))Re(hA(x))
Im(hA(x))

0
-5
-10

-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10x

Figure 8: Eq. (4.17) に至るまでの計算過程. X = [−10, 10], N = 200, K = 1/∆x ≈ 10 として, プログラム (6) を用いて計算. 図中 fA (k), gA (k), hA (k) は, 解析解 (4.14) を表す. 解は Eq. (3.80)Eq. (3.82)

44

4 プログラム

4.3

DFT の計算速度

4.3.1

DFT の高速化

4.3 DFT の計算速度

さて, DFT のコード (1) を見ると, 1 回の DFT ごとに Descrepter 等を作成しています.
複数回 DFT を実行する場合, これは一回だけで良いはずなのでその分計算時間が掛かってしまいそうです. 以下のコードを考えて比較してみましょう. 基本的にはコード (1) と同じですが, サブルーチン dft の中を取り出して main の中に移しています. (2) の中に, mainの中の重要な部分だけを取り出して記述しました. また, 並列計算の影響もなくすために

1 スレッドだけで比較しています.
Listing 2: 離散フーリエ変換の点数 Np = 218 として 10 × 100 回実行した時の計算時間を測定するプログラム
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

call mkl_set_num_threads(1) ! Number of threads for MKL parallel computation! DFT preparation with double-complex type of size N array
Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,Np)
Status = DftiSetValue(hand,DFTI_FORWARD_SCALE,1d0)
Status = DftiSetValue(hand,DFTI_BACKWARD_SCALE,1d0/dble(Np))
Status = DftiCommitDescriptor(hand)!-------------------------f0=ff1=ftime_fast=0e0
time_simp=0e0
do j=1,100
write(6,*)j
time_fast0=0e0
time_simp0=0e0
time_fast1=0e0
time_simp1=0e0
call cpu_time(time_simp0)
do i=1,10
call dft(Np,f0,"forward")
call dft(Np,f0,"backward")
enddo
call cpu_time(time_simp1)
time_simp = time_simp + time_simp1-time_simp0
call cpu_time(time_fast0)
do i=1,10
Status = DftiComputeForward(hand,f1)
Status = DftiComputeBackward(hand,f1)
enddo
call cpu_time(time_fast1)
time_fast = time_fast + time_fast1-time_fast0
enddo!----------------------------! Free DFT settings
Status = DftiFreeDescriptor(hand)
write(6,*)"time simple [CPUsec] ",time_simp
write(6,*)"time fast [CPUsec] ",time_fast

実行結果はListing 3: 計算速度比較結果
1
2

time simple [CPUsec] 5.71431065
time fast [CPUsec] 4.00213337

45

4 プログラム

Figure 9:

4.3 DFT の計算速度同一計算機上でのベンチマーク結果.

コンパイラ：ifort (IFORT) 16.0.2 20160204,

MKL:/opt/intel/compilers and libraries 2016.2.181, 計算に用いたスレッド数は 1.

となりました. 可読性のために使用した dft 関数を使わなければ, 約 0.7 倍の計算時間となり高速化されます. 可読性を重視するか速度を重視するかは状況次第ですので, ここでは結果を出すにとどめておきます.

4.3.2

DFT の分点数に対する計算速度

MKL のルーチン（MKL は結局 ﬀtw と呼ばれるルーチンを利用していますが…）は,任意の数に対してフーリエ変換を実行することが可能です. つまり, 高速フーリエ変換だ
けサポートしているのではなく, 分点数が奇数でも素数でも使用可能です. ですが, 分点数
がちょうど 2, 3, 5, 7, 11 の乗数に等しいときに自動的に高速フーリエ変換が実行されます.
計算速度を調べたのが Fig. 9 です. 横軸に分点数, 縦軸に計算時間を示しています.

4.3.3

並列計算時の比較

さて, MKL は並列計算をサポートしています. 離散フーリエ変換も例外ではなく, 並列に CPU を使用することによって高速化することができます. 並列計算の準備は以下の一行プログラム内で指定するだけで良く, あとは勝手に実行されます.
Listing 4: MKL ルーチン実行時に使用するスレッド数の設定
1

call mkl_set_num_threads(Nthreads)46

4 プログラム

4.3 DFT の計算速度

Figure 10: 同一計算機上でのベンチマーク結果. 計算時間の CPU 個数 (8 コア 16 スレッド) と配列の要素
数依存性.

プログラムが始まったらすぐに指定すると良いでしょう. もしくは, MKL ルーチンの実行
直前が良いでしょう.

Fig. 10 に, 並列計算の速度, 配列のサイズと実行時間の関係を示します.
計算 1 回 とは（順方向→逆方向フーリエ変換）の 1 セットを 1 回 としています.

8 コア 16 スレッドの環境で実行したので, 図の場合では物理的な CPU の個数 8 個で計算速度はおおよそ打ち止めになっています.

47

REFERENCE

REFERENCE

Reference[1] インテル (R) マス・カーネル・ライブラリーリファレンスマニュアル (2006 年)[2] W. H. Press, S. A. Teukolsky, W. T. Vetterling and Brian P. Flannery, (丹慶勝市,
奥村晴彦, 佐藤俊郎, 小林誠 訳), 『ニューメリカルレシピ・イン・シー 日本語版―C
言語による数値計算のレシピ』, 技術評論社 (1993)
英語かつ過去版で良ければ, 公式から無料で閲覧することができます (http://

numerical.recipes/oldverswitcher.html).

48

A フーリエ変換→逆変換で元に戻ることの証明

A

フーリエ変換→逆変換で元に戻ることの証明

A.1

連続の場合

∫ ∞

dk̃f (k̃)ei2πk̃x
[∫ ∞
]
∫−∞
∞
′
′ −i2π k̃x′=dk̃dx f (x )eei2πk̃x
−∞
∫−∞
∫ ∞
∞
′
′
′
=
dx f (x )dk̃ei2πk̃(x−x )−∞∫−∞∞dx′ f (x′ )δ(x − x′ )
=

(A.1)

= f (x)(A.5)

f (x) =

(A.2)
(A.3)
(A.4)

−∞

A.2

離散の場合

N
−1
∑

=

=

=

=

m=0
N
−1
∑
m=0
N
−1
∑
n=0
N
−1
∑
n=0
N
−1
∑

e−i2πkn′ (xm −xoﬀ ) · f (xm )
e−i2πkn′ (xm −xoﬀ ) ·N−1∑

f (kn )ei2πkn (xm −xoﬀ )(A.6a)(A.6b)

n=0

f (kn )f (kn )

N
−1
∑
m=0
N
−1
∑

ei2π·(kn −kn′ )·(xm −xoﬀ )′

ei2π·(n−n )∆k·m∆x(A.6c)(A.6d)

m=0

f (kn )N δn,n′(A.6e)

n=0= N · f (kn′ )(A.6f)となります.Eq. (4.7) と, 等比級数の知識であるN−1(a = 0)
 N,
∑
i2πna
e
= ei2πaN − 1, (a =
̸ 0)n=0
ei2πa − 1

(A.7)

を用いて,
N
−1
∑

′

ei2πm(n−n )/N = N δn,n′

m=0

49

(A.8)

B DFT の具体例で考えた無理数倍についてを用いています. ここで, δi,j はクロネッカーのデルタで 1, (i = j)
δi,j =
 0, (i ̸= j)

(A.9)

です.

B

DFT の具体例で考えた無理数倍について
項 3.1.6 で考えた, 離散フーリエ変換の区間が波の周期の無理数倍を仮定したことについてコメントします. 今考えたい計算は以下の通りです.

f (kn ) =
=

N
−1
∑
m=0
N
−1
∑

f (xm )e−i2πnm/N
cos(2π3 · m

m=0
N
−1 [
∑

(B.1a)

X −i2πnm/N)eN](3X−n)(3X+n)
1ei2πm N + e−i2πm N
2 m=0( i2π(3X−n))e−11e−i2π(3X+n) − 1=+
2 ei2π(3X−n)/N − 1 e−i2π(3X+n)/N − 1( i2π3N ∆x)e−11e−i2π3N ∆x − 1=+2 e−i2π(kn −3)∆x − 1 e−i2π(kn +3)∆x − 1=(B.1b)(B.1c)...?

(B.1d)(B.1e)となります. Eq. (B.1c) から Eq. (B.1d) に移る際に (3X ± n)/N はゼロではないという条件を使用しています. 右辺の分母に注目すると, kn = ±3 の時にどちらかの分母がゼロとなり発散, 計算が不能になります. つまり, Eq. (B.1c) ではどう頑張っても N が有限ならば値も有限ですが, Eq. (B.1d) では無限があり得るといっています. 明らかにおかしいですので, 本来 Eq. (B.1d) への変形は行うことができません.
つまりこの条件が言っているのはどんなに偶然でも kn = n∆k が 3 になった場合にはEq. (B.1e) を使用することはできないということです. すなわち, f (x) の周期 Xp の整数倍となるときは適当な整数 n′ を用いて

kn = n∆k = n

1
n′
=
X
Xp

(B.2)

が満たされるときは, Eq. (B.1e) は使えない, X について X が Xp の整数倍ではないときに

X=

n
Xp
n′

(B.3)

は使えないことを意味します. 以上をまとめると, Eq. (B.1e) が使用できるのは, 区間

x = [0, X] の X は波の波数の無理数倍じゃなきゃダメです, ということです.

50

B DFT の具体例で考えた無理数倍について
実際に離散フーリエ変換を行う上で無理数倍の区間指定はほとんど行いません. そのため, 具体例として挙げるのは適切ではないと判断し, 項 3.1.6 ではこれ以上の式変形を行っておりません.

51

C 任意区間の離散フーリエ変換プログラム

C

任意区間の離散フーリエ変換プログラム
Listing 5: Intel MKL を用いた Fortran90 のプログラム例 (位置を範囲変更した一般の場合)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

program main
implicit none
integer::mdouble precision::dx,Xrange,xa,xb
double precision,allocatable::x(:)
integer::ndouble precision::dk,Krange
double precision,allocatable::k(:)
integer::Npcomplex(kind(0d0)),allocatable::f(:)
complex(kind(0d0)),external::func
Np = 200 ! Number of DFT points
xa = -exp(1d0)xb = exp(1d0)
Xrange = xb-xa ! x=[xa, xa + Xrange]! Initialize and allocation for DFT
allocate(x(0:Np-1),k(0:Np-1),f(0:Np-1))
x = 0d0
k = 0d0
f = dcmplx(0d0,0d0)! Define position x
dx = Xrange / Np
call dft_abscissa_shift(Np, x, dx, xa) ! x=[xa, xa + Xrange]! Calculate k abscissa
Krange = 1d0/dx
dk = Krange / Np
call dft_abscissa(Np, k, dk) ! k=[0, Krange]! Calculate function f at x=x(m)do m = 0,Np-1
f(m) = func(x(m))
enddo! Write down f(x)call dft_write("before_fx.d",Np,x,f)! Forward dft
call dft(Np,f,"forward")! Write down f(k)call dft_writefk_shift("fk.d",Np,k,f,xa) ! output k=[0,Krange]
call dft_writefk_center0_shift("fk_center0.d",Np,k,f,xa) ! output k=[-Krange/2,Krange/2]! Backward dft
call dft(Np,f,"backward")! Write down f(x)call dft_write("after_fx.d",Np,x,f)
stop
end program main
complex(kind(0d0)) function func(x)implicit none
double precision,intent(in)::x
double precision::pi=dacos(-1d0)
func=dcmplx(dcos(2d0*pi*3d0*x),0d0)
returnend function func!---------------------------------subroutine dft_abscissa(N, w, dw)
implicit none
integer,intent(in)::N
double precision,intent(in)::dw
double precision,intent(out)::w(0:N-1)

52

C 任意区間の離散フーリエ変換プログラム

75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

!
!w=[0, W]! w(0) = 0! w(N-1) = (N-1)*dw = W-dw!
integer::j
do j = 0,N-1
w(j) = j*dw
enddo
returnend subroutine dft_abscissa
subroutine dft_abscissa_shift(N, w, dw, wa)
implicit none
integer,intent(in)::N
double precision,intent(in)::dw,wa
double precision,intent(out)::w(0:N-1)!
!w=[wa, wa + W]
! w(0) = wa
! w(N-1) = (N-1)*dw = W - dw + wa!
integer::j
do j = 0,N-1
w(j) = j*dw + wa
enddo
returnend subroutine dft_abscissa_shift
include "/opt/mkl/include/mkl_dfti.f90"
subroutine dft(N,f,FB)
use MKL_DFTI
implicit none
integer,intent(in)::N
complex(kind(0d0)),intent(inout)::f(0:N-1)
character(*),intent(in)::FB!
!sikinote! author : sikino! date : 2022/06/11!
!n: Number of DFT points.
!f(i): value of data at i
!FB: "forward" -> Forward DFT! "backward" -> Backward DFT
integer::Status
TYPE(DFTI_DESCRIPTOR),POINTER::hand
Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,N)
Status = DftiSetValue(hand,DFTI_FORWARD_SCALE,1d0)
Status = DftiSetValue(hand,DFTI_BACKWARD_SCALE,1d0/dble(N))
Status = DftiCommitDescriptor(hand)if(trim(FB)=="forward")then
Status = DftiComputeForward(hand,f)elseif(trim(FB)=="backward")then
Status = DftiComputeBackward(hand,f)else
write(6,*)"DFT string different"
stop
endif
Status = DftiFreeDescriptor(hand)
returnend subroutine dft
subroutine dft_write(fname,N,w,f)
implicit none
character(*),intent(in)::fname
integer,intent(in)::N
double precision,intent(in)::w(0:N-1)
complex(kind(0d0)),intent(in)::f(0:N-1)! w = [wa,W]! w(0) = wa
! --> output range w=[wa,wa+W]
integer::j53

C 任意区間の離散フーリエ変換プログラム

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225

open(21,file=trim(fname))
do j = 0,N-1
write(21,’(3e26.16e3)’)w(j),dble(f(j)),dimag(f(j))
enddoclose(21)returnend subroutine dft_write
subroutine dft_writefk_shift(fname,N,k,f,xa)
implicit none
character(*),intent(in)::fname
integer,intent(in)::N
double precision,intent(in)::k(0:N-1),xa
complex(kind(0d0)),intent(in)::f(0:N-1)! k = [0,K]! k(0) = 0! --> output range k=[0,K]
integer::nncomplex(kind(0d0))::tmp
double precision,parameter::pi=dacos(-1d0)
open(21,file=trim(fname))
do nn = 0,N-1
tmp = exp(dcmplx(0d0,2d0*pi*k(nn)*xa))*f(nn)write(21,’(3e26.16e3)’)k(nn),dble(tmp),dimag(tmp)
enddoclose(21)returnend subroutine dft_writefk_shiftsubroutine dft_writefk_center0_shift(fname,N,k,f,xa)
implicit none
character(*),intent(in)::fname
integer,intent(in)::N
double precision,intent(in)::k(0:N-1),xa
complex(kind(0d0)),intent(in)::f(0:N-1)! k = [0,K]! k(0) = 0! --> output range k=[-K/2,K/2]
integer::nn,Nhalf
double precision::dk,Krange
complex(kind(0d0))::tmp
double precision,parameter::pi=dacos(-1d0)
if(mod(N,2)==0)then
Nhalf = N/2
else
Nhalf = (N-1)/2
endif
dk = abs(k(1)-k(0))
Krange = N*dk
open(21,file=trim(fname))
do nn = Nhalf,N-1
tmp = exp(dcmplx(0d0,2d0*pi*k(nn)*xa))*f(nn)write(21,’(3e26.16e3)’)k(nn)-Krange,dble(tmp),dimag(tmp)
enddo
do nn = 0,Nhalf-1
tmp = exp(dcmplx(0d0,2d0*pi*k(nn)*xa))*f(nn)write(21,’(3e26.16e3)’)k(nn),dble(tmp),dimag(tmp)
enddoclose(21)returnend subroutine dft_writefk_center0_shift54

D 任意区間の離散フーリエ変換を利用する畳み込みプログラム

D

任意区間の離散フーリエ変換を利用する畳み込みプログラム
いくつかのサブルーチンは Program (5) と共通なので, 省略しています.
Listing 6: Intel MKL を用いた畳み込みのプログラム例 (位置を範囲変更した一般の場合)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

program main
implicit none
integer::mdouble precision::dx,Xrange,xa,xb
double precision,allocatable::x(:)
integer::ndouble precision::dk,Krange
double precision,allocatable::k(:)
integer::Npcomplex(kind(0d0)),allocatable::f(:),g(:),h(:)
complex(kind(0d0)),external::func,gfunc
double precision,parameter::pi=dacos(-1d0)
Np = 200 ! Number of DFT points
xa = -10d0
xb = 10d0
Xrange = xb - xa ! x=[xa,xa+Xrange]! Initialize and allocation for DFT
allocate(x(0:Np-1),k(0:Np-1),f(0:Np-1))
x = 0d0
k = 0d0
f = dcmplx(0d0,0d0)
allocate(g(0:Np-1),h(0:Np-1))
g = dcmplx(0d0,0d0)
h = dcmplx(0d0,0d0)! Define position x
dx = Xrange / Np
call dft_abscissa_shift(Np, x, dx, xa) ! x=[xa, xa + Xrange]! Calculate k abscissa
Krange = 1d0/dx
dk = Krange / Np
call dft_abscissa(Np, k, dk) ! k=[0, Krange]! Calculate function f at x=x(m)do m = 0,Np-1
f(m) = func(x(m))g(m) = gfunc(x(m))
enddo! Write down f(x)call dft_write("before_fx.d",Np,x,f)call dft_write("before_gx.d",Np,x,g)! Forward dft
call dft(Np,f,"forward")
call dft(Np,g,"forward")! Multiply f(k) and g(k) and add the phase factor in k-space
do n = 0,Np-1h(n) = exp(dcmplx(0d0,-2*pi*k(n)*xa))*f(n)*g(n)enddo! Write down f(k), g(k), h(k)call dft_writefk_center0_shift("fk.d",Np,k,f,xa) ! output k=[-Krange/2,Krange/2]
call dft_writefk_center0_shift("gk.d",Np,k,g,xa) ! output k=[-Krange/2,Krange/2]
call dft_writefk_center0_shift("hk.d",Np,k,h,xa) ! output k=[-Krange/2,Krange/2]! Backward dft
call dft(Np,h,"backward")! Write down h(x)call dft_write("after_hx.d",Np,x,h)stop
end program main
complex(kind(0d0)) function func(x)implicit none55D 任意区間の離散フーリエ変換を利用する畳み込みプログラム

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

double precision,intent(in)::x
func=dcmplx(x**5*exp(-x**2),0d0)
returnend function func
complex(kind(0d0)) function gfunc(x)implicit none
double precision,intent(in)::x
gfunc=dcmplx(exp(-4*x**2),0d0)
returnend function gfunc

56