フーリエ解析#
フーリエ解析(Fourier analysis):信号を異なる周波数の正弦波(sin)の重ね合わせとして表現すること
例えば画像を周波数の表現にして、周波数をフーリエ解析で分解したりする。
フーリエ級数#
幅\(T\)の区間\([-T/2, T/2]\)における連続関数\(f(t)\)はフーリエ級数に展開できる
フーリエ係数は
である。 ここで
である。これを 基本周波数 といい、フーリエ級数は関数\(f(t)\)を\(\omega_o\)の整数倍の周波数の正弦波の重ね合わせとして表現する。
\(a_0/2\)を 直流成分 と呼び、\(\omega_o\)の\(k\)倍の周波数の振動を第\(k\) 高調波 と呼ぶ。以下では\(t\)を「時刻」とみなし、\(f(t)\)を「信号」と呼ぶ。
複素数の指数関数#
複素数#
複素数とは \(z=x+i y\) のように 実(数)部 \(x\) に 虚(数)部 \(y\) の \(i\) 倍を足したもの
\(i\) は 虚数単位 と呼ばれ、 \(i^2=-1\) となる数と約束したもの。
\(x\)軸を実部 \(x=\operatorname{Re} z\)、\(y\)軸を 虚部 \(y=\operatorname{Im} z\)にとった平面を 複素平面 と呼ぶ
Show code cell source
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(4, 3))
ax.arrow(-0.5, 0, 1.7, 0, head_width=0.05, head_length=0.1, fc='black', ec='black', length_includes_head=True)
ax.arrow(0, -0.5, 0, 1.5, head_width=0.05, head_length=0.1, fc='black', ec='black', length_includes_head=True)
ax.text(1.3, 0, "Re", fontsize=14, color='black')
ax.text(0, 1.1, "Im", fontsize=14, color='black')
ax.set_aspect('equal', adjustable='datalim')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.set(xticks=[], yticks=[])
c = 1 + 0.5j # 複素数
complex_points = [c, c.conjugate()] # conjugate: 共役複素数
colors = ['red', 'blue']
for i, (point, color) in enumerate(zip(complex_points, colors)):
x, y = point.real, point.imag
ax.plot([0, x], [0, y], color=color, linestyle='--', linewidth=0.8)
ax.scatter([x], [y], color=color)
z = r"$z$" if i == 0 else r"$\bar{z}$"
ax.text(x + 0.05, y + 0.05, f'{z} = ${point}$', fontsize=10, color=color)
plt.show()
共役複素数#
複素数 \(z=x+i y\) に対して, \(\bar{z}=x-i y\) をその 共役(きょうやく)複素数 といい、 \(z, \bar{z}\) は 互いに複素共役である という。 これらは複素平面上で \(x\) 軸に関して対称な位置にある。 この定義から
となる。また重要な関係として
が成り立つ。
オイラーの式#
指数部が虚数の指数関数 \(e^{i \theta}\) を次の複素数と定義する。
これは オイラーの式 と呼ばれる式で、複素平面の単位円上の実軸から角度\(\theta\)の点を表す。
が成り立つ。
証明
定義より次のように変形できる
両者が等しいことは、三角関数の加法定理
よりわかる
が成り立つ。
証明
オイラーの式 およびその式の \(\theta\) を \(-\theta\) に置き換えた式
を使うと、
より
が得られる。また
より
が得られる。
フーリエ級数の複素表示#
を使うと、フーリエ級数は以下のように表すことができる。
ただし
である
また、\(C_k\)の各式は次のように変形できる
この\(C_k\)を (複素)フーリエ係数 と呼ぶ。
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(4, 3))
# 軸を描画
ax.arrow(-1.2, 0, 2.4, 0, head_width=0.05, head_length=0.1, fc='black', ec='black', length_includes_head=True)
ax.arrow(0, -1.2, 0, 2.4, head_width=0.05, head_length=0.1, fc='black', ec='black', length_includes_head=True)
ax.text(1.3, 0, "Re", fontsize=14, color='black')
ax.text(0, 1.3, "Im", fontsize=14, color='black')
ax.set_aspect('equal', adjustable='datalim')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.set(xticks=[], yticks=[])
# 単位円
theta = np.linspace(0, 2*np.pi, 100)
x = np.cos(theta)
y = np.sin(theta)
plt.plot(x, y, 'k-', linewidth=1) # 黒い実線で単位円を描画
# curve
t = 0.3*np.pi
theta = np.linspace(0, t, 100)
x = np.cos(theta) * 0.2
y = np.sin(theta) * 0.2
ax.plot(x, y, 'k--', linewidth=1) # 黒い実線で単位円を描画
# text
m = len(x) // 2
ax.text(x[m], y[m], r"$\theta$")
# line
x = np.cos(t)
y = np.sin(t)
ax.plot([0, x], [0, y], color="blue", linestyle='-.', linewidth=0.8)
ax.scatter([x], [y], color="blue")
plt.show()
フーリエ変換#
フーリエ級数
において
とおき、
と書くと、
となる。周期\(T\to\infty\)とすると、\(\Delta \omega \to 0\)となり、積分として表すことができる。
を信号\(f(t)\)の フーリエ変換 と呼び、
を 逆フーリエ変換 と呼ぶ。
Note
信号\(f(t)\)(逆フーリエ変換)は周波数\(\omega\)の正弦波\(e^{i \omega t}\)を負の数を含めたすべての実数\(\omega\)に対して重ね合わせている
フーリエ変換 \(F(\omega)\) は信号\(f(t)\)に\(e^{i \omega t}\)の共役複素数\(e^{-i \omega t}\)を掛けて、すべての時刻\(t\)で積分した形になっている
フーリエ変換
逆フーリエ変換
逆フーリエ変換は信号\(f(t)\)をあらゆる周波数(すべての実数の\(\omega\))の振動の重ね合わせで表す。
\(F(\omega)\)は周波数\(\omega\)の成分\(e^{i \omega t}\)の大きさを表し、\(f(t)\)の スペクトル (フランス語: spectre、英語: spectrum)と呼ばれる。(光をプリズムに通して各色のスペクトルに分けるイメージ)
\(|\omega|\)が大きい成分は 高周波成分 、小さい成分は 低周波成分 と呼ばれる。とくに\(\omega=0\)の成分は定数であり、 直流成分 と呼ばれる。
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
# サンプルデータを生成する
N = 1000
d = 0.0001 # サンプリング周期
t = np.arange(0, N * d, d) # 時間
def f(t):
freq1, freq2 = 100, 500 # 周波数
s1, s2 = 1.5, 2 # スペクトル F(ω)
return s1 * np.sin(2 * np.pi * t * freq1) + s2 * np.sin(2 * np.pi * t * freq2)
y = f(t)
F = np.fft.fft(y) # フーリエ変換
freq = np.fft.fftfreq(N, d=d) # 周波数スケール
# 振幅スペクトル Amplitude の取得
# 振幅スペクトルは信号をフーリエ変換した結果の絶対値をとったもの
Amp = np.abs(F)
Amp = Amp / (N / 2) # 正規化
# 結果をプロット
fig, axes = plt.subplots(nrows=2)
axes[0].plot(t, y)
axes[0].set(title=r'Original Function $f(t)$', xlabel=r"$t$", ylabel=r"f(t)")
# 左右対称なので正の値だけ N//2 で取り出す
axes[1].stem(freq[:N//2], Amp[:N//2], 'r', markerfmt=" ", basefmt=" ")
axes[1].set(title=f'Fourier Transform $F(\omega)$', xlabel=r'Frequency $ω$', ylabel=r'Magnitude $F(ω)$')
fig.tight_layout()
fig.show()
例:矩形窓#
次の関数は幅\(W\)の 矩形窓(rectangular window) と呼ばれる
このフーリエ変換は
となる。ただし、
と定義する。
Show code cell source
# 参考:numpyの高速フーリエ変換 fft を使う場合
# 自分で離散的なデータポイントからフーリエ変換するのは離散フーリエ変換を学んでからやる
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(4, 3))
# もとの関数 f(t)
w = 1.5 # width
def f(t):
return (1/2*w) * 1 * ((-w <= t) & (t <= w))
t = np.linspace(-3, 3, 100)
y = f(t)
ax.plot(t, f(t), label="f(t)", color="black", alpha=0.8)
# フーリエ変換 F(omega)
fourier = np.fft.fft(y)
inv_fourier = np.fft.ifft(fourier).real # 実部のみ
ax.plot(t, inv_fourier, label=r"$\hat{f}(f)$", color="red", linestyle="--", alpha=0.7)
ax.set(xlabel="t", ylabel="f(t)")
ax.legend()
plt.show()
次の関数\(w(t)\) を幅\(\sigma\)の ガウス窓(Gaussian window) という。
このフーリエ変換は
ここで \(z=t+i \sigma^2 \omega\)と変数変換すると
となる
Show code cell source
from scipy import signal
import matplotlib.pyplot as plt
window = signal.windows.gaussian(51, std=7)
fig, ax = plt.subplots(figsize=(4, 2))
ax.plot(window)
ax.set(title=r"Gaussian window ($\sigma$=7)", ylabel="Amplitude", xlabel="Sample")
fig.show()
Show code cell source
import numpy as np
from scipy import signal
from scipy.fft import fft, fftshift
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(4, 2.4))
window = signal.windows.gaussian(51, std=7)
A = fft(window, 2048) / (len(window)/2.0)
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
ax.plot(freq, response)
ax.plot(window)
ax.set(title=r"Frequency response of the Gaussian window ($\sigma$=7)",
ylabel="Normalized magnitude [dB]", xlabel="Normalized frequency [cycles per sample]",
xlim=(-0.5, 0.5), ylim=(-120, 0))
fig.show()
畳み込み積分#
たたみこみ積分
次の積分を信号 \(f(t), g(t)\) の たたみこみ積分 (convolution) または 合成積 と呼ぶ
(左辺の積分における\(f\)と\(g\)の変数の和が\(s + t-s = t\)となっている)
畳み込み積分は分配法則、交換法則、結合法則が成り立つ
分配法則: \(f(t) *(a g(t)+b h(t))=a f(t) * g(t)+b f(t) * h(t)\)
交換法則:\(f(t) * g(t)=g(t) * f(t)\)
結合法則:\((f(t) * g(t)) * h(t)=f(t) *(g(t) * h(t))\)
定理(信号の合成積はフーリエ変換の積)
信号 \(f(t), g(t)\) のフーリエ変換をそれぞれ \(F(\omega), G(\omega)\) とするとき、 \(f(t) * g(t)\) のフーリエ変換は \(F(\omega) G(\omega)\) である
証明
\(f(t) * g(t)\)のフーリエ変換は次のようになる
(積分の順序を入れ換え、 \(t^{\prime}=t-s\) と変数変換し、 \(t=t^{\prime}+s\) を代入した)
一般に、ある集合に要素間の演算が定義されているものを 代数系 と呼ぶ。そして2つの異なる集合での演算が同じ規則に従うとき、2つの代数系は 同型 であるという。 関数の和\(+\)と畳み込み積分\(*\)に関する演算は、そのフーリエ変換の和\(+\)と積\(\times\)に関する演算と同じであり、両方の代数系が同型である。
フィルター#
信号 \(f(t)\) の値を時刻 \(t\) をはさむ幅 \(2 W\) の区間で平均し、これを \(\tilde{f}(t)\) と定義する。
矩形窓やガウス窓\(w(t)\)を使うと、次のように書き直すことができる
これにより、関数\(w(t)\)を様々な関数に変えれば信号\(f(t)\)の様々な変換ができる。このような操作を関数\(w(t)\)による フィルター という。
\(w(t), f(t)\)のフーリエ変換をそれぞれ\(W(\omega), F(\omega)\)とおけば、畳み込み積分の定理により\(\tilde{f}(t) = w(t)*f(t)\)のフーリエ変換は
と書くことができる。
もし\(|\omega|\)が小さい部分で\(|W(\omega)|\)が大きく、\(|\omega|\)が大きい部分で\(|W(\omega)|\)が小さくなる場合、\(w(t)\)によるフィルターは低周波成分を増幅し、高周波成分を減衰させる。このようなフィルターを ローパスフィルター ( 低域フィルター )と呼ぶ。その逆のフィルターを ハイパスフィルター ( 高域フィルター )と呼ぶ。
その中間として、特定の\(|\omega|\)の付近を増幅し、それ以外を減衰させるフィルターは 帯域フィルター ( バンドパスフィルター )と呼ぶ。
Show code cell source
import numpy as np
import matplotlib.pyplot as plt
N = 1000
t = np.linspace(0, 10, N)
def f(t):
return np.sin(t) + 0.5 * np.sin(10 * t)
signal = f(t)
fig, ax = plt.subplots(figsize=(4, 3))
ax.plot(t, signal, label="f(t)", color="black", alpha=0.8)
ax.set(xlabel="t", ylabel="f(t)")
from scipy.signal.windows import gaussian
from scipy.signal import fftconvolve
def gaussian_smoothing(signal, window_size, std_dev):
# ガウス窓を作成
gauss_window = gaussian(window_size, std=std_dev)
# gauss_window /= gauss_window.sum() # 窓の正規化
# 畳み込み
return fftconvolve(signal, gauss_window, mode='same') # 同じ長さを保つ
smoothed_signal = gaussian_smoothing(signal, window_size=N//2, std_dev=0.25)
ax.plot(t, smoothed_signal, label=r"$\tilde{f}(t)$", color="steelblue", alpha=0.8)
ax.legend()
fig.show()
小まとめ
信号\(f(t)\)をフーリエ変換
して、個々の\(\omega\)に対応するスペクトル(周波数\(\omega\)の影響度の強さ)\(F(\omega)\)を得ることができ、またそれにより\(f(t)\)に含まれる周波数ごとの正弦波\(e^{i\omega t}\)も分解できる。
そして窓関数\(w(t)\)による畳み込み
をすることで、特定の周波数をカットするなどの変換を行って信号を変換することができる。
一つの身近な具体例は画像に対するガウシアンフィルタ(ガウシアン平滑化・ガウシアンブラー)によるぼかし
参考#
金谷健一. (2003). これなら分かる応用数学教室: 最小二乗法からウェーブレットまで.