MDで数式を扱う練習

WARNING

教養課程レベルの行列/微積分学しか知らないので、測度論みたいな厳密な話はあんまりしない。

INFO

Fourier変換の大体の話はInternation Tables for Crystallography Vol.B 1.3に書いてある。

フーリエ変換 (Fourier Transform)

Fourier Transform

$$\begin{array}{rl}F(k)& ={\int }_{-\infty }^{\infty }f(x)e^{-2\pi ikx}dx\\ f(x)& ={\int }_{-\infty }^{\infty }F(k)e^{2\pi ikx}dk\end{array}$$

3Blue1BrownのBut what is the Fourier Transform? A visual introduction.が幾何的イメージでフーリエ変換を表現していてとてもよい。複素平面上の大きさ１の円周上に変換したい関数を巻き付けている。

fourier_plot.jl

using Plots, FFTW

sin2pi(x) = sinpi(2x)
cos2pi(x) = cospi(2x)

xs = [exp(-x^2) * (sin2pi(7x) + 0.7sin2pi(16x) + 1.5cos2pi(32x)) for x in range(-2, 2, 512)]

## 4/512 = 0.0078125 is the pixel size
## 1/(2*0.0078125) = 64 is the Nyquist frequency

xs_fft_pw = xs |> rfft .|> abs2

## Plot both the original and the power spectrum in two subplots

plot(
    plot(xs, title="Signal", legend=false, xticks=(0:64:512, -2:0.5:2)),
    plot(xs_fft_pw, title="Power spectrum", legend=false, xticks=(0:32:256, 0:8:64)),
    layout=(2, 1), size=(600, 300)
)

fourier_plot2.jl

# ...Continue from the previous code
fourier_curve(k) = xs .* cis.(-2pi * k .* range(-2, 2, length=N)) ##　フーリエ変換の積分の中

data(k) = fourier_curve(k) |> x -> (real(x), imag(x))
point(k) = fourier_curve(k) / N*4 |> sum |> x -> (real(x), imag(x))　## ただ足すとデカくなりすぎるから適当に割ってる

pl1 = plot(data(15), title="k = 15", aspect_ratio=:equal, legend=false, size=(600, 600))
scatter!(point(15))
pl2 = plot(data(16), label="",title="k = 16", aspect_ratio=:equal, legend=false, size=(600, 600))
scatter!(point(16), label="")
plot(pl1, pl2, layout=(2, 1), size=(600, 600))

フーリエ変換で実空間の関数を逆空間に変換した図。複素数はそのままだとプロットし辛いのでよくパワースペクトラム（絶対値の二乗）をプロットする（CTFの図とか）。位相の情報が失われていることに注意！ フーリエ変換の積分の中身をそのままプロットした図。オレンジの点が係数になってる（大きさはみやすさのために適当）。k=15ではほぼ０だけど、k=16になると下の方にズレて値を持つようになる。これは位相の情報も残ってる。

これらが幾何的にわかりやすい説明。正規直交基底で係数を抜き出すイメージもあるけど、省略。

フーリエ変換されるものは実数値関数である必要はなく、$f(x)\in {\mathbb{C}}^n$でもよい。実際一般に$F(k)\in {\mathbb{C}}^n$なので、双対性より$f(x)\in {\mathbb{C}}^n$となってもおかしくない。 （だいたいのセンサがパワーを計測してるので）画像とか音声を扱っていると勘違いしがち。

畳み込みの定理　(Convolution Theorem)

Convolution Theorem

$$\begin{array}{rl}\{f\ast g\}(x)& ={\int }_{-\infty }^{\infty }f(t)g(x-t)dt\\ \mathcal{F}[f\ast g]& =\mathcal{F}[f]\mathcal{F}[g]\end{array}$$

証明

$$\begin{array}{rl}\mathcal{F}[f\ast g]& ={\int }_{-\infty }^{\infty }({\int }_{-\infty }^{\infty }f(t)g(x-t)dt)e^{-2\pi ikx}dx\\ & ={\int }_{-\infty }^{\infty }f(t)({\int }_{-\infty }^{\infty }g(x-t)e^{-2\pi ikx}dx)dt\\ & ={\int }_{-\infty }^{\infty }f(t)({\int }_{-\infty }^{\infty }g(x)e^{-2\pi ik(x+t)}dx)dt\\ & ={\int }_{-\infty }^{\infty }f(t)({\int }_{-\infty }^{\infty }g(x)e^{-2\pi ikx}dx)e^{-2\pi ikt}dt\\ & =\mathcal{F}[f]\mathcal{F}[g]\end{array}$$

（積分の順序って交換してええんかな）

フーリエ変換でいちばん大事な性質。実空間での積は逆空間の畳み込みになり、逆空間での積は実空間の畳み込みになる。この性質はめちゃくちゃ多用する。

さっきの図は実空間でexp(-x^2)を掛けていたけど、exp(-64x^2)に変えたもの。パワースペクトラムがより滲んでいるのがわかる。

イメージとしては、$f$のあらゆる場所$x$に$f(x)$で重みをつけた$g$を置いて足し合わせたもの。例えば格子点に$f$を畳み込むと、格子点周りに$f$が出現する。散乱理論とかグリーン関数が出てくるものでこのイメージを使う。（ちゃんと散乱理論をやったことはない）この性質は次の定理でつかう。

サンプリング定理 (Nyquist–Shannon sampling theorem)

Nyquist–Shannon sampling theorem

関数$f(x)$のフーリエ変換$F(k)$が$F(k)=0(|k|\ge f_{\mathrm{nyq}})$となる周波数$f_{\mathrm{nyq}}$を持つとき、$f(x)$を間隔$\mathrm{\Delta }=1/(2f_{\mathrm{nyq}})$でサンプリングしたデータから復元できる。

補題（くし関数）

$$\begin{array}{rl}{\text{Ш}}_{\mathrm{\Delta }}(t)& =\sum _{n=-\infty }^{\infty }\delta (t-n\mathrm{\Delta })\\ \mathcal{F}[{\text{Ш}}_{\mathrm{\Delta }}(t)]& =\frac{1}{\mathrm{\Delta }}{\text{Ш}}_{1/\mathrm{\Delta }}(k)\\ & =\frac{1}{\mathrm{\Delta }}\sum _{n=-\infty }^{\infty }\delta (k-n/\mathrm{\Delta })\end{array}$$

実空間で$\mathrm{\Delta }$間隔でデルタ関数が周期的にあると、そののフーリエ変換は逆空間で$1/\mathrm{\Delta }$間隔でデルタ関数が周期的にある。

解説

前述の畳み込みの定理と、くし関数の性質を使う。

定理の$F(k)$の条件を満たす関数$f(x)$を考える。$f(x)$を$\mathrm{\Delta }$間隔でサンプリング（実空間上での積）すると、逆空間上では$F(k)$を$1/\mathrm{\Delta }$間隔に置いた畳み込みになる。

このとき、逆空間上では$\pm f_{\mathrm{nyq}}$の範囲に$F(k)$があるので、実空間でサンプリングを$1/(2f_{\mathrm{nyq}})$間隔で行うと、逆空間では$2f_{\mathrm{nyq}}$間隔で$F(k)$を置くので、周波数が重ならないで情報が保存される。

前提の$F(k)=0(|k|\ge f_{\mathrm{nyq}})$が成り立たないと、$F(k)$が重なってしまうので情報が混じってしまう。これをエイリアシングという。

具体例でいうと、音声を44.1kHzでサンプリングしておけば、人間の聴覚範囲（20kHz）までの音声を復元できる。単粒子解析では目標分解能$1/f_{\mathrm{res}}$Åに対応してピクセルサイズが$1/(2f_{\mathrm{res}})$Åになるようにサンプリングする。

デジタルカメラのセンサーの前に、よく光学ローパスフィルター（OLPF）が入っている。これは、センサーのナイキスト周波数を超える周波数の光をカットするためのもので、モアレを防ぐために必要。でも最近の高画素化で現実的に問題となる高周波信号はないとか、後処理でローパスをかけたりできるからOLPFがない機種もある。そういう光学素子は厚みに依存して軸外の光学特性が悪くなりそうだし納得はできる。Sigmaの高画素機fp LはOLPFがない。

電子顕微鏡の場合、ピクセルサイズを決定したらそれ以上の高周波信号はノイズになるので、対物絞りでローパスフィルターをかける場合がある。しかし、CCPEMのメーリス曰く、対物絞りの金属板のコンディションによりいらん収差を生むことがあるので、入れないほうがいいかも。（電子回折以外で使ったことない）最近だとFalcon4とかが超解像度モードを搭載していて、8KモードでPhysical Nyquistを超えた分解能が出る場合もあるから、対物絞りは使わないが吉。

光学顕微鏡の構造化照明法は、サンプルに照明する光のパターンを高周波の波にすることで（実空間での積）、逆空間の高分解能な情報を得る手法。あんまり詳しくは知らないけどこの原理を使って超解像を実現している。電顕でもレーザー位相板の技術を使ってこれ出来ないかなーとか妄想してる。

スーパーヘテロダイン受信機の原理もこれと同じだと思うけど、詳しくは知らない。

長くなったので次回に続く。