9. フーリエ変換

あらゆる信号には、時間領域でも周波数領域でも解析できるスペクトルがあります。あるいは、空間領域や空間周波数領域での測定も可能です。フーリエ変換は、信号を時間領域から周波数領域へ変換する数学的操作です。

私たちは物質の光学応答を解析するためにフーリエ変換を使います。物質の時間依存非線形応答を測定したとき、振動するシグネチャが現れることがあります。これらの振動の周波数は、系内のコヒーレンスや結合に対応します。

また、フーリエ変換赤外 (FTIR) 分光器という装置もよく使われます。この装置は、サンプルから来た光の干渉パターン (インターフェログラム) を測定し、フーリエ変換を使ってその干渉パターンをスペクトルに変換します。

フーリエ変換をきちんと扱おうとすれば何週間もかかります。ここでは、Julia でフーリエ変換を計算し、その結果を解釈するために必要な基本だけを扱います。

2 つの周波数

60 Hz の正弦波と、その基本周波数の半分の振幅を持つ第 1 高調波 120 Hz の正弦波を考えましょう。下の図のように時間領域の信号は複雑に見えますが、この信号のフーリエ変換を取ると、60 Hz と 120 Hz に明瞭な 2 つのピークが現れます。この簡単な例は、複雑な信号をその構成周波数に分解するというフーリエ変換の威力を示しています。光と物質の相互作用から生じるはるかに複雑な信号の解析にも、この手法を使います。

これをコードで実装しましょう。まず必要なパッケージをインポートします。可視化用の GLMakie とフーリエ変換用の FFTW です。 (FFTW.jl は、C で書かれた FFTW ライブラリ (Fastest Fourier Transform in the West) のラッパーにすぎません。) そして信号関数を定義します。

using FFTW
using GLMakie

signal(t, f1, f2) = sin(2π * f1 * t) + 0.5 * sin(2π * f2 * t)

次にサンプリング周波数を決め、時間ベクトルを生成します。サンプリング周波数とは、1 秒あたりに信号をサンプルする回数のことです。ここではサンプリング周波数の逆数をサンプル間の時間間隔として定義します。サンプリング周波数の半分はナイキスト周波数と呼ばれ、正確に表現できる最大の周波数です。信号を最高周波数の 2 倍より低いレートでサンプリングすると、エイリアシングが発生し、信号が実際よりも低い周波数に見えてしまいます。ぜひ低いレートでサンプリングしてみて、何が起こるか観察してみてください。

fs = 1000 # sample rate in Hz
t_max = 1 # max time in seconds
t = 0:1/fs:t_max

f1 = 60 # frequency in Hz
f2 = 120 # frequency in Hz
y = signal.(t, f1, f2)

次のコードで信号のフーリエ変換を計算できます。

Y = abs.(fftshift(fft(y)))
X = fftshift(fftfreq(length(t), fs))

信号とそのフーリエ変換のプロットがうまくいったら、コード中の 3 つの関数 — fft、fftshift、fftfreq — がそれぞれ何をしているかを考えてみましょう。

fft 関数は、信号の離散フーリエ変換 (DFT) を計算量効率の良い方法で計算します。結果は複素数の配列なので、振幅スペクトルを得るために絶対値を取っています。

fftfreq 関数は、ウィンドウサイズ n = length(t) (時間ベクトルの長さ) とサンプリング間隔 fs をもとに、フーリエ変換に対応する周波数を計算します。正の周波数を先に、続けて負の周波数を返します。フーリエ変換は 0 を中心に対称で、区間 $(-\infty, \infty)$ で定義されます。したがって、フーリエ変換の後半は前半の鏡像にすぎません。 fftshift 関数は、ゼロ周波数成分をスペクトルの中央に移動させ、可視化を容易にします。 fftfreq の例は、Julia REPL で ?fftfreq と入力するか、VS Code で関数名にマウスを乗せれば見ることができます。同様に fft の結果もスペクトルの中央に揃えるためにシフトが必要です。

fftshift の効果は、1 から 10 までの数列に対して呼び出すとよくわかります。

julia> fftshift(1:10)
10-element Vector{Int64}:
  6
  7
  8
  9
 10
  1
  2
  3
  4
  5

配列の前半と後半を入れ替えているだけだとわかります。信号のフーリエ変換を fftshift を使った場合と使わない場合でプロットしてみて、見た目を比較してみてください。

問題

デルタ関数 $\delta(x)$ は、次の性質を持つ興味深い関数です。

\delta(x) = \begin{cases} 0, & x \ne 0, \\ \infty, & x = 0, \end{cases}

そして、

\int_{-\infty}^{\infty} \delta(x)\, dx = 1.

ディラックのデルタ関数は実際には関数ではなく分布 (distribution) で、分布は積分の性質、あるいはテスト関数に対する作用によって定義されます。 $x = 0$ で連続な関数 $f(x)$ に対して、デルタ関数は次の性質を持ちます。

\int_{-\infty}^{\infty} f(x)\, \delta(x)\, dx = f(0),

$x = a$ で連続なら、次が成り立ちます。

\int_{-\infty}^{\infty} \delta(x - a)\, f(x)\, dx = \int_{-\infty}^{\infty} \delta(a - x)\, f(x)\, dx = f(a).

そして、 $\delta(x - a)$ のフーリエ変換は次のようになります。

\int_{-\infty}^{\infty} \delta(x - a)\, e^{-i x \xi}\, dx = e^{-i a \xi}.

ディラックのデルタの定義を使って、次を計算してみましょう。

\frac{1}{2} \int_{-\infty}^{\infty} [\delta(x + a) + \delta(x - a)]\, e^{-i x \xi}\, dx

時定数と振動周波数を使って減衰振動子の式を書き、それをコードで実装しましょう。続いてフーリエ変換を計算します。両方の曲線をプロットしてください。時定数を変えたとき、時間領域の曲線とその変換のどちらにも何が起こるでしょうか?

以下は矩形波とそのフーリエ変換を実装したコードです。まずコードを実行せず、読むだけで 以下の問いに答えてください。

矩形波の周波数はいくつでしょうか?
フーリエ変換にはいくつの高調波が現れますか? それらの周波数は?
周波数分解能はどれくらいでしょうか?
与えられたコードのナイキスト周波数はいくつでしょうか? 何番目の高調波の後に現れますか? つまり、正確に表現できる最大の周波数は?

ここで コードを実行 し、以下の問いに答えてください。

高調波の数を増やすとフーリエ変換はどう変わりますか? 観察した結果を答えてください。次に点数を減らしてみる (つまりナイキスト周波数を下げる) と、どうなるでしょうか?

タンパク質中の光合成集光アンテナは電子コヒーレンスを持つことが知られていますが、生体環境のノイズの多さから測定が極めて難しいことでも有名です。この問題で扱う論文は、結合したクロモフォア間のエキシトンと分子振動の間の長寿命コヒーレンスを測定しています。以下の問題は、論文の Figure 2d の生データを使います。生データと基本的なコードはチュートリアルのリポジトリにあります。解析を進めて、振動の寿命と周波数を抽出してください。

Zhu, et al. Quantum Phase Synchronization via Exciton-Vibrational Energy Dissipation Sustains Long-Lived Coherence in Photosynthetic Antennas. Nat Commun 2024, 15 (1), 3171. https://doi.org/10.1038/s41467-024-47560-6.

単一指数減衰関数を使ってデータをフィッティングし、エネルギー緩和時間を抽出してください (フィッティングの章の NonlinearCurveFitProblem と solve を使います)。データとフィット結果をプロットしてください。
データからフィット結果を引き算し、残差のフーリエ変換を計算してください。フーリエ変換をプロットし、振動の周波数を特定してください。この周波数の物理的な意味は何でしょうか?