宿題 1 のグラフ描き

画像処理とフーリエ変換
宿題 1 のグラフ描き
桂田祐史
http://nalab.mind.meiji.ac.jp/~mk/fourier/
2014 年 10 月 24 日
問 11 の関数 f , g の Fourier 級数展開そのものは授業中に説明した通り、
(
)
π
4 cos x cos 3x cos 5x
f (x) = −
+
+
+ ··· ,
2 π
12
32
52
(
)
4 sin x sin 3x sin 5x
+
+
+ ··· .
g(x) =
π
1
3
5
f は連続で区分的に C 1 級、g は区分的に C 1 級 (ゆえに 2 乗可積分) であるから、どのような
収束をするかは、10/17 の授業で説明した 3 つの定理で説明が出来る。
こうすれば Fourier 級数の部分和 fn , gn のグラフが描ける
f[n_,x_]:=Pi/2-4/Pi Sum[Cos[k x]/k^2,{k,1,n,2}]
gf10=Plot[f[10, x], {x, -10, 10}]
g[n_,x_]:=4/Pi Sum[Sin[k x]/k,{k,1,n,2}]
gg20=Plot[g[20,x],{x,-10,10}]
n を変えて gn の変化の様子を見たければ、
Manipulate[Plot[g[n, x], {x, -Pi, Pi}, PlotPoints -> 200], {n, 1, 100, 1}]
グラフィックスを保存するには
Export["graphf10.pdf",gf10]
Export["graphg20.pdf",gg20]
g の不連続点 x = 0, ±π の近くでは、ぴょこっと高い山、それと対称な低い谷があり、n が
増加すると横幅は小さくなるが、高さはほとんど変わらない (低くなってくれない) ことが見
て取れる。
おまけ: hn のグラフを描いてみる
h[n_,x_]:=4/Pi Sum[Cos[k x],{k,1,n,2}]
Manipulate[Plot[h[n,x], {x,-7,7}, PlotRange -> Full, AspectRatio->Automatic,
PlotPoints -> 100], {n, 2, 20, 2}]
, AspectRatio->Automatic は入れた方が良いかどうか判らないが。
1
http://nalab.mind.meiji.ac.jp/~mk/lecture/fourier-2014/toi1.pdf
1
もっと色々 Mathematica にさせてみる
f や g のグラフを描かせたり、Fourier 係数を求めて、部分和の計算をさせたりしてみよう。
f (x) = |x| (−π < x < π) ということで、積分の計算をするだけならば
f[x_]:=Abs[x]
使うだけで、十分であるが、周期関数らしいグラフを描くには次のようにすると良い。
f[x_]:=Abs[Mod[x,2Pi,-Pi]]
gf=Plot[f[x],{x,-10,10}]
とすると良い。Fourier 係数 ak =
1
π
∫π
−π
f (x) cos kx dx の計算は、例えば
Integrate[f[x]Cos[k x],{x,-Pi,Pi}]/Pi
これで
2(−1+cos kπ+kπ sin kπ)
k2 π
が得られる。k が整数と教えてないので、結果が複雑。
Simplify[Integrate[f[x]Cos[k x],{x,-Pi,Pi}]/Pi, Assumptions->Element[k,Integers]]
2(−1+(−1)k )
とすると、
k2 π
が得られる。実はこれにも問題がある。k = 0 のとき正しくない式で
ある！
！ここら辺は Mathmatica は変な割り切り方をしている (マニュアルによると、“不定積
分において，Integrate は，ほとんどすべてのパラメータの値に対して正しい結果を求める．”
— 少数の例外は無視？)。注意が必要である。k = 0 の場合については別扱いで
Integrate[f[x],{x,-Pi,Pi}]/Pi,
とすると良い。π という結果が得られる。
∫π
bk = π1 −π f (x) sin kx dx については
Integrate[f[x] Sin[k x], {x, -Pi, Pi}]
とすると 0 が返る。
f の Fourier 係数を a[] という関数で計算させて、f の n 項までの部分和 fn を計算するよ
うにする。
ほぼ全自動
Clear[a,f]
f[x_]:=Abs[Mod[x,2Pi,-Pi]]
a[k_]=Simplify[Integrate[f[x]Cos[k x],{x,-Pi,Pi}]/Pi, Assumptions->Element[k,Integers]]
a[0]=Integrate[f[x],{x,-Pi,Pi}]/Pi
f[n_,x_]:=a[0]/2+Sum[a[k]Cos[k x],{k,1,n}]
gf10=Plot[f[10,x],{x,-10,10}]
(実は Simplify[] はしなくてもあまり効率に影響がない。f[n,x] は Sum[] を使っているの
で、遅延評価 := をしないとうまく動かない。一方 a[k ] の方を遅延評価するのはお勧めでき
ない。ここら辺は難しい。最初に Fourier 係数を全部計算して、その結果をリストに保存して
使うのが良いかもしれない。)
2
x = 0 のところで fn (x) − f (x) の値を調べると、図 1 のように減少している (1/n に比例し
ていることが分かる)。
list=Table[f[2^k,0.0]-f[0.0],{k,1,10}]
gerror=ListLogPlot[list]
Export["error.pdf",gerror]
n = 210 のとき誤差は 0.000621699 ≒ 0.6 × 10−3
1
0.50
0.10
0.05
10-2
10-3
2
4
6
8
10
図 1: fn (0) − f (0) の変化 (n = 2k , k = 1, 2, · · · , 10)
Clear[g,bg]
g[x_]=Sign[Mod[x,2Pi,-Pi]]
bg[k_]=Simplify[Integrate[g[x]Sin[k x],{x,-Pi,Pi}]/Pi,
Assumptions->Element[k,Integers]]
g[n_,x_]:=Sum[bg[k]Sin[k x],{k,1,n}]
ggandg10=Plot[{g[x],g[10,x]},{x,-10,10}]
1.0
0.5
-10
5
-5
-0.5
-1.0
図 2: g と g の 10 項までの部分和 g10 のグラフ
3
10
h について
h(x) = 2
∑
(−1)n δ(x − nπ).
n∈Z
ここで δ は Dirac のデルタ関数 (以下では単にデルタ関数という) を表す。デルタ関数は
実は関数ではなく、超関数というものであり、本来はグラフを描いたりすることも出来ない
ものであるが (x ̸= 0 で 0, x = 0 で ∞ という「説明」もあるが、∞ は図に描けない)、実は
Fourier 級数展開は可能であり、その部分和は以下に見るように普通の関数なので、グラフを
描くことが出来る。
デルタ関数は、任意の連続関数 φ に対して
∫ ∞
δ(x)φ(x) dx = φ(0)
−∞
を満たすという条件や、
∫
δ(x) = 0 (x ̸= 0),
∞
δ(x) dx = 1
−∞
で特徴づけられる。物理学では、大きさのない 1 点に質量や電荷などが集中しているという理
想化を考えるが (質量の場合は質点、電荷の場合は点電荷と呼ぶ)、デルタ関数は、そのように
1 点に 1 単位のものが集中している場合の、密度を表す関数と考えられる (直観的には、原点
で無限大。積分すると総量になって 1)。
h の定義には、デルタ関数を平行移動した δ(x − nπ) が現れるが、δ(x − a) は
∫ ∞
δ(x − a)φ(x) dx = φ(a)
−∞
という性質を持つ。
さて、h の Fourier 級数展開の計算であるが、f や g のように [−π, π] で積分するとして、
どうすれば良いだろうか？
(♯1 )
h(x) = 2 (δ(x) − δ(x − π)) ,
(♯2 )
h(x) = 2 (δ(x) − δ(x + π)) ,
(♯3 )
h(x) = 2 (δ(x) − δ(x − π) − δ(x + π)) ,
このどれが正しいのだろうか？周期 2π の周期関数とは R/2πZ で定義された関数のことと考
えると、π と −π は同じ点を表すことになるので、(♯1 ), (♯2 ) のどちらも同じことを表してい
て正しそうであり、(♯3 ) は 1 つのものを重複してカウントしていて間違っていそうである。
あるいは積分範囲を [−π, π] でなく、少し右にずらした [−π + ε, π + ε] を考えてみると (ε
は小さい正の数)、(♯1 ) を採用することになり、逆に [−π − ε, π − ε] と左にずらすと、(♯2 ) を
採用することになる。どちらで計算しても Fourier 係数は同じになる。
これから
∫
)
2
2(
1 π
h(x) cos kx dx = (cos (k · 0) − cos kπ) =
1 − (−1)k ,
ak =
π −π
π
π
∫ π
1
2
bk =
h(x) sin kx dx = (sin (k · 0) − sin kπ) = 0.
π −π
π
4
ゆえに h の Fourier 級数展開は
h(x) =
4
(cos x + cos 3x + cos 5x + · · · ) .
π
この級数は普通の意味では収束しない、発散級数であるが、超関数の意味では収束している。
また部分和は普通の関数になっていて、グラフも描ける。
h[n_,x_]:=4/Pi Sum[Cos[k x],{k,1,n,2}]
Manipulate[Plot[h[n,x], {x,-7,7}, PlotRange -> Full, AspectRatio->Automatic,
PlotPoints -> 100], {n, 2, 20, 2}]
工学系の本には「デルタ関数列」というものが良く出て来る (数学の解析の本で、軟化作用
素を定義する際の関数列と同じものと考えて良い)。これは n → ∞ の極限で δ となる関数列
{φn } である。
∫
∞
−∞
φn (x) dx = 1
という性質を満たしつつ、n が大きくなるにつれ、x = 0 の近くに集中していくものになって
いる。例えば工学の本で良く見かけるのは
φn (x) =
sin nx
.
πx
このグラフも描いてみて、hn のグラフと見比べてみると良い。
なお、ここまで来ると推察できると思うが、h は g の (超関数の意味での) 導関数である。
実際、Heaviside の階段関数
{
1 (x ≥ 0)
H(x) =
0 (x < 0)
の導関数は δ になることが知られているので、
g ′ (x) = 2δ(x) − 2δ(x − π) = h(x) ([−π + ε, π + ε] での等式).
そして、これまでと同様に、h の Fourier 級数展開は、g の Fourier 級数展開を項別微分した
ものになっている。
5

Download Report