juliaで再学習フーリエ変換（その0）

理工系界隈でフーリエ変換とよく耳にしますが……

フーリエ変換という単語を聞きき、理工系の方でこの単語を知らない方は少ないと思います。

僕自信の経験談として、修士のときまで、偏微分方程式を専攻していましたが……
結局、フーリエ変換は、偏微分方程式を解く道具としか、当時考えておらず……
音声処理や画像処理に使われているという事実だけ、知っている程度で、どのように応用されているか知りませんでした。
社会人になり、エンジニアになったときに「高速フーリエ変換」という単語を知りました。
当時、「高速」ってなに？
「フーリエ変換に早いも遅いもあるの？」と単語に純粋に疑問を感じたことだけあります。
今は、科学計算系の会社のエンジニア？になって、フーリエ変換を学部で習う計算問題から、数学科が道具として使うときや音声処理や画像処理などについて、まとめてみようと思い立ちこの記事を書きました。
言語は、話題性、早い、実行環境がJupyterなど、手軽で話題性があるという点で、Julia v1.1.0で書いていきます。

とりあえず、フーリエ変換級数展開

フーリエ変換と聞いて、まず思い出すのがフーリエ級数です。

区間 $[-\pi, \pi]$ 上の関数 $f$ を以下のように展開する

f = \displaystyle {\frac {a_{0}}{2} +\sum _{k=1}^{\infty }(a_{k}\cos kx+b_{k}\sin kx)} \frac{a_0}{2} + \sum_{k=1}^\infty (a_k\cos kx+b_k\sin kx) \\
\begin{eqnarray}
    \begin{cases}
        \displaystyle a_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f\left(t \right) \cos nt\,dt \ \left(n = 0,1,2,3,\cdots \right) \notag \\
        \displaystyle b_n = \frac{1}{\pi} \int_{-\pi}^{\pi} f\left(t \right) \sin nt\,dt \ \left(n = 0,1,2,3,\cdots \right)
    \end{cases}
\end{eqnarray}

ことをフーリエ変換展開といいます。

関数にいろいろ制約は付きますが、とりあえず、初回なので、固いことは、後に回しましょう。

簡単な例

\begin{eqnarray}
	f(x) =
	\begin{cases}
		x & (0 \leqq x \leqq 1) \notag \\
		0 & (-1 \leqq x \leqq 0)
	\end{cases}
\end{eqnarray}

function f(x)
	if x >= 0. && x <= 1.
		val = x
	elseif x >= -1. && x <= 0.
		val = 0.
	else
		val = 0.
	end
	return val
end

フーリエ級数展開すると

f(x) = \frac{1}{4} - \frac{2}{\pi^2} \left( \cos{\pi x} + \frac{\cos{3\pi x}}{3^2} + \cdots + \frac{\cos{(2n-1)\pi x}}{(2n-1)^2} + \cdots \right) + \frac{1}{\pi} \left( \sin{\pi x} - \frac{\sin{2\pi x}}{2} + \cdots + (-1)^{n-1} \frac{\sin{n\pi x}}{n} + \cdots \right)

function f_fourier_seri(x, n)
	val = 1/4
	for i in 1:n
		val -= 2/(π^2) * (cos.((2*i-1)*π*x)/(2*i-1)^2)
		val += 1/π * ((-1)^(i-1) * sin.((i*π*x))/i)
	end
	return val
end

描画すると

x = -1.5:0.1:1.5
y0=[f(xx) for xx in x]
anim = @animate for i=1:1:100
    y1=[f_fourier_seri(xx, i) for xx in x]
    plot(x,y0, color="red", alpha = 0.8)
    plot!(x,y1, color="darkcyan", alpha = 0.4) 
end
gif(anim, "output_gif/fourier_anim_fps60.gif", fps = 60)

Gist