フーリエ解析

フーリエ解析#

フーリエ解析（Fourier analysis）：信号を異なる周波数の正弦波（sin）の重ね合わせとして表現すること

例えば画像を周波数の表現にして、周波数をフーリエ解析で分解したりする。

フーリエ級数#

幅Tの区間[−T/2,T/2]における連続関数f(t)はフーリエ級数に展開できる

f(t)=a02+∑k=1∞(akcos⁡kωot+bksin⁡kωot),−T2≤t≤T2

フーリエ係数は

ak=2T∫−T/2T/2f(t)cos⁡kωot dt,bk=2T∫−T/2T/2f(t)sin⁡kωot dt

である。 ここで

ωo=2πT

である。これを 基本周波数 といい、フーリエ級数は関数f(t)をωoの整数倍の周波数の正弦波の重ね合わせとして表現する。

a0/2を 直流成分 と呼び、ωoのk倍の周波数の振動を第k 高調波 と呼ぶ。以下ではtを「時刻」とみなし、f(t)を「信号」と呼ぶ。

Tip

基本周波数 ω0=2π/Tは、幅T、区間[−2/T,2/T]を1周期とする振動の周波数である。

例えばsin⁡(x)は1周期が2πなのでω0=1となる

複素数の指数関数#

複素数#

複素数とは z=x+iy のように 実（数）部 x に 虚（数）部 y の i 倍を足したもの

i は 虚数単位 と呼ばれ、 i2=−1 となる数と約束したもの。

x軸を実部 x=Rez、y軸を 虚部 y=Imzにとった平面を 複素平面 と呼ぶ

Hide 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()
../../_images/04c0e0d48a3b3947862799efb52199d5ea79623c6064b59696dea22ec0e50051.png

共役複素数#

複素数 z=x+iy に対して, z¯=x−iy をその 共役（きょうやく）複素数 といい、 z,z¯ は 互いに複素共役である という。 これらは複素平面上で x 軸に関して対称な位置にある。 この定義から

x=z+z¯2,y=z−z¯2i

となる。また重要な関係として

zz¯=(x+iy)(x−iy)=x2+y2=|z|2

が成り立つ。

オイラーの式#

指数部が虚数の指数関数 eiθ を次の複素数と定義する。

eiθ=cos⁡θ+isin⁡θ

これは オイラーの式 と呼ばれる式で、複素平面の単位円上の実軸から角度θの点を表す。

定理
eiθeiϕ=ei(θ+ϕ)

が成り立つ。

証明

定義より次のように変形できる

eiθeiϕ=(cos⁡θ+isin⁡θ)(cos⁡ϕ+isin⁡ϕ)=(cos⁡θcos⁡ϕ−sin⁡θsin⁡ϕ)+i(cos⁡θsin⁡ϕ+sin⁡θcos⁡ϕ)ei(θ+ϕ)=cos⁡(θ+ϕ)+isin⁡(θ+ϕ)

両者が等しいことは、三角関数の加法定理

cos⁡(θ+ϕ)=cos⁡θcos⁡ϕ−sin⁡θsin⁡ϕsin⁡(θ+ϕ)=cos⁡θsin⁡ϕ+sin⁡θcos⁡ϕ

よりわかる

定理
cos⁡θ=eiθ+e−iθ2,sin⁡θ=eiθ−e−iθ2i

が成り立つ。

証明

オイラーの式 およびその式の θ を −θ に置き換えた式

eiθ=cos⁡θ+isin⁡θ,e−iθ=cos⁡θ−isin⁡θ

を使うと、

eiθ+e−iθ=cos⁡θ+isin⁡θ+cos⁡θ−isin⁡θ=2cos⁡θ

より

cos⁡θ=eiθ+e−iθ2

が得られる。また

eiθ−e−iθ=cos⁡θ+isin⁡θ−cos⁡θ+isin⁡θ=2isin⁡θ

より

sin⁡θ=eiθ−e−iθ2i

が得られる。

フーリエ級数の複素表示#

cos⁡θ=eiθ+e−iθ2,sin⁡θ=eiθ−e−iθ2i

を使うと、フーリエ級数は以下のように表すことができる。

f(t)=a02+∑k=1∞(akeikωot+e−ikωot2+bkeikωot−e−ikωot2i)=a02+12∑k=1∞((ak−ibk)eikωot+(ak+ibk)e−ikωot)=∑k=−∞∞Ckeikωot

ただし

Ck={(ak−ibk)/2k>0a0/2k=0(a−k+ib−k)/2k<0

である

また、Ckの各式は次のように変形できる

ak−ibk2=2T∫−T/2T/2f(t)cos⁡kωot−isin⁡kωot2 dt=1T∫−T/2T/2f(t)e−ikωot dta02=1T∫−T/2T/2f(t)dta−k+ib−k2=2T∫−T/2T/2f(t)cos⁡(−kωot)+isin⁡(−kωot)2 dt=1T∫−T/2T/2f(t)(cos⁡kωot−isin⁡kωot)dt=1T∫−T/2T/2f(t)e−ikωot dt
f(t)=∑k=−∞∞Ckeikωot,Ck=1T∫−T/2T/2f(t)e−ikωot dt

このCkを （複素）フーリエ係数 と呼ぶ。

Hide 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()
../../_images/3ac9a136235eef71f0516135228216d7d80bcbbbc37ddc69f827571591f58d95.png

フーリエ変換#

フーリエ級数

f(t)=∑k=−∞∞Ckeikωot,Ck=1T∫−T/2T/2f(t)e−ikωot dt

において

ωk=kωo,k=0,±1,±2,…

とおき、

Δω=ωk+1−ωk=ωo=2πT,Ck=F(ωk)T

と書くと、

f(t)=12π∑k=−∞∞F(ωk)eiωktΔω,F(ωk)=∫−T/2T/2f(t)e−iωkt dt

となる。周期T→∞とすると、Δω→0となり、積分として表すことができる。

F(ω)=∫−∞∞f(t)e−iωt dt

を信号f(t)の フーリエ変換 と呼び、

f(t)=12π∫−∞∞F(ω)eiωt dω

を 逆フーリエ変換 と呼ぶ。

Note

  • 信号f(t)（逆フーリエ変換）は周波数ωの正弦波eiωtを負の数を含めたすべての実数ωに対して重ね合わせている

  • フーリエ変換 F(ω) は信号f(t)にeiωtの共役複素数e−iωtを掛けて、すべての時刻tで積分した形になっている

フーリエ変換

F(ω)=∫−∞∞f(t)e−iωt dt

逆フーリエ変換

f(t)=12π∫−∞∞F(ω)eiωt dω

逆フーリエ変換は信号f(t)をあらゆる周波数（すべての実数のω）の振動の重ね合わせで表す。

F(ω)は周波数ωの成分eiωtの大きさを表し、f(t)の スペクトル （フランス語: spectre、英語: spectrum）と呼ばれる。（光をプリズムに通して各色のスペクトルに分けるイメージ）

|ω|が大きい成分は 高周波成分 、小さい成分は 低周波成分 と呼ばれる。とくにω=0の成分は定数であり、 直流成分 と呼ばれる。

図：フーリエ変換のイメージ。F(ω)は一般に複素数なのでグラフに書けない → 絶対値をとる

Hide 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'Absolute of Fourier Transform $|F(\omega)|$', xlabel=r'Frequency $ω$', ylabel=r'Magnitude $|F(ω)|$')

fig.tight_layout()
fig.show()
../../_images/8d35437b0f81d22db1fd74a98f10a007b331fb02c5041c5c8c4df5a08abca68d.png

例：矩形窓#

次の関数は幅Wの 矩形窓（rectangular window） と呼ばれる

w(t)={1/2W−W≤t≤W0 それ以外 

このフーリエ変換は

W(ω)=12W∫−WWe−iωt dt=12W∫−WW(cos⁡ωt−isin⁡ωt)dt(オイラーの式よりe−iωt=cos⁡ωt−isin⁡ωt)=12W(∫−WWcos⁡ωt dt−∫−WWisin⁡ωt dt)(∵積分の線形性)=12W∫−WWcos⁡ωt dt(∵sin⁡は奇関数すなわちsin⁡(−θ)=−sin⁡(θ)であり、−Wから+Wまで積分したら0になるため消える)=1W∫0Wcos⁡ωt dt(∵∫−WWcos⁡ωt dt=2∫0Wcos⁡ωt dt)=1W[sin⁡ωtω]0W(∵∫cos⁡(ωt)dt=sin⁡(ωt)ω)=1W(sin⁡ωWω−sin⁡0ω)=sin⁡WωWω(∵sin⁡0=0)=sincWπω

となる。ただし、

sincx:=sin⁡πxπx

と定義する。

Hide 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()
../../_images/e0749b41f0acbdbf87fedc38b75cb73b404d389a5898a5502a1ea22b681e3517.png

次の関数w(t) を幅σの ガウス窓（Gaussian window） という。

w(t)=12πσe−t2/2σ2

このフーリエ変換は

W(ω)=∫−∞∞w(t)e−iωt dt=12πσ∫−∞∞e−t2/2σ2e−iωt dt=12πσ∫−∞∞e−(t22σ2+iωt) dt(指数法則anam=an+m)=12πσ∫−∞∞e−(t+iσ2ω)2/2σ2−σ2ω2/2 dt(平方完成したり色々いじるとこうなる)=(12πσ∫−∞∞e−(t+iσ2ω)2/2σ2 dt)e−σ2ω2/2

ここで z=t+iσ2ωと変数変換すると

W(ω)=(12πσ∫−∞∞e−z2/2σ2 dz)e−σ2ω2/2=e−σ2ω2/2=2πσ(12πσ−1e−ω2/2σ−2)

となる

Hide 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()
../../_images/f4e461a241b6ad84b835617f393b3c5c5ecf604d9bb9311c9491cbc5ae1b72c6.png
Hide 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()
../../_images/23b1609dcf0186b8605a0e594b8ddfd99dfa7f7ffb2c2d568659c7681e9e39f6.png

（出典：gaussian — SciPy v1.14.1 Manual ）

畳み込み積分#

たたみこみ積分

次の積分を信号 f(t),g(t) の たたみこみ積分 （convolution） または 合成積 と呼ぶ

f(t)∗g(t)=∫−∞∞f(s)g(t−s)ds

（左辺の積分におけるfとgの変数の和がs+t−s=tとなっている）

畳み込み積分は分配法則、交換法則、結合法則が成り立つ

  • 分配法則： f(t)∗(ag(t)+bh(t))=af(t)∗g(t)+bf(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(ω),G(ω) とするとき、 f(t)∗g(t) のフーリエ変換は F(ω)G(ω) である

証明

f(t)∗g(t)のフーリエ変換は次のようになる

∫−∞∞f(t)∗g(t)e−iωt dt=∫−∞∞(∫−∞∞f(s)g(t−s)ds)e−iωt dt=∫−∞∞f(s)(∫−∞∞g(t−s)e−iωt dt)ds=∫−∞∞f(s)(∫−∞∞g(t′)e−iω(t′+s)dt′)ds=∫−∞∞f(s)(∫−∞∞g(t′)e−iωt′dt′)e−iωs ds=(∫−∞∞f(s)e−iωs ds)(∫−∞∞g(t′)e−iωt′dt′)=F(ω)G(ω)

（積分の順序を入れ換え、 t′=t−s と変数変換し、 t=t′+s を代入した）

一般に、ある集合に要素間の演算が定義されているものを 代数系 と呼ぶ。そして2つの異なる集合での演算が同じ規則に従うとき、2つの代数系は 同型 であるという。 関数の和+と畳み込み積分∗に関する演算は、そのフーリエ変換の和+と積×に関する演算と同じであり、両方の代数系が同型である。

フィルター#

信号 f(t) の値を時刻 t をはさむ幅 2W の区間で平均し、これを f~(t) と定義する。

f~(t):=12W∫−WWf(t−s)ds

矩形窓やガウス窓w(t)を使うと、次のように書き直すことができる

f~(t)=∫−∞∞w(s)f(t−s)ds=w(t)∗f(t)

これにより、関数w(t)を様々な関数に変えれば信号f(t)の様々な変換ができる。このような操作を関数w(t)による フィルター という。

w(t),f(t)のフーリエ変換をそれぞれW(ω),F(ω)とおけば、畳み込み積分の定理によりf~(t)=w(t)∗f(t)のフーリエ変換は

F~(ω)=W(ω)F(ω)

と書くことができる。

もし|ω|が小さい部分で|W(ω)|が大きく、|ω|が大きい部分で|W(ω)|が小さくなる場合、w(t)によるフィルターは低周波成分を増幅し、高周波成分を減衰させる。このようなフィルターを ローパスフィルター （ 低域フィルター ）と呼ぶ。その逆のフィルターを ハイパスフィルター （ 高域フィルター ）と呼ぶ。

その中間として、特定の|ω|の付近を増幅し、それ以外を減衰させるフィルターは 帯域フィルター （ バンドパスフィルター ）と呼ぶ。

例

fftconvolve — SciPy v1.14.1 Manual

ガウシアンブラーで画像を平滑化（ぼかす）例がある

Hide 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=100)
ax.plot(t, smoothed_signal, label=r"$\tilde{f}(t)$", color="steelblue", alpha=0.8)

# scipy
from scipy.ndimage import gaussian_filter1d
smoothed_signal = gaussian_filter1d(signal, 100)
ax.plot(t, smoothed_signal, label=r"$\tilde{f}(t)$ (scipy gaussian_filter1d)", color="orange", alpha=0.8)

ax.legend()
fig.show()
../../_images/4dc6dc747004cc62dd8578347e88dec2c460a1fd2077e08ea7fb51d08bbdc353.png

小まとめ

信号f(t)をフーリエ変換

F(ω)=∫−∞∞f(t)e−iωt dt

して、個々のωに対応するスペクトル（周波数ωの影響度の強さ）F(ω)を得ることができ、またそれによりf(t)に含まれる周波数ごとの正弦波eiωtも分解できる。

そして窓関数w(t)による畳み込み

f~(t)=w(t)∗f(t)=∫−∞∞w(s)f(t−s)ds

をすることで、特定の周波数をカットするなどの変換を行って信号を変換することができる。

一つの身近な具体例は画像に対するガウシアンフィルタ（ガウシアン平滑化・ガウシアンブラー）によるぼかし

パワースペクトル#

パーセバル（・プランシュレル）の式

信号 f(t),g(t) のフーリエ変換をそれぞれ F(ω),G(ω) とする。

このとき、以下の パーセバル（・プランシュレル）の式 が成り立つ

∫−∞∞f(t)g(t)― dt=12π∫−∞∞F(ω)G(ω)― dω∫−∞∞|f(t)|2 dt=12π∫−∞∞|F(ω)|2 dω
証明
∫−∞∞f(t)g(t)― dt=12π∫−∞∞F(ω)G(ω)― dω

については、

∫−∞∞f(t)g(t)―dt=∫−∞∞f(t)(12π∫−∞∞G(ω)―e−iωt dω)dt=12π∫−∞∞(∫−∞∞f(t)e−iωt dt)G(ω)―dω=12π∫−∞∞F(ω)G(ω)―dω

g(t)=f(t)とすると

∫−∞∞|f(t)|2 dt=12π∫−∞∞|F(ω)|2 dω

が得られる

信号f(t)が時間tとともに変動する振動を表すとき、|f(t)|2はその振動の単位時間あたりの エネルギー とみなすことができる。

パーセバルの第2式

∫−∞∞|f(t)|2 dt=12π∫−∞∞|F(ω)|2 dω

はエネルギーがすべてのωにわたる|F(ω)|2の積分で表されることを意味する。

これを パワースペクトル と呼び、

P(ω)=|F(ω)|2

は周波数ωの振動成分のもつエネルギーであると解釈できる。

フーリエ変換F(ω)は一般に複素数だが、その絶対値は実数になりグラフに描くことができる。

フーリエ変換は共役F(−ω)=F(ω)―であり|F(−ω)|=|F(ω)―|だから、パワースペクトルはグラフにすると原点を中心に左右対称になる。

Hide code cell source 
# うまくいかなかった
import numpy as np
import matplotlib.pyplot as plt

sampling_rate = 30
duration = 10
N = int(sampling_rate * duration)

t = np.linspace(0, duration, N, endpoint=False)
frequencies = [30, 60, 90, 150, 300]
amplitudes = [1.0, 0.8, 0.6, 0.4, 0.2]
signal = sum(
    amp * np.sin(2 * np.pi * freq * t) for freq, amp in zip(frequencies, amplitudes)
)

# Perform the Fourier Transform
fft_result = np.fft.fft(signal)
frequencies = np.fft.fftfreq(len(signal), 1 / sampling_rate)
power_spectrum = np.abs(fft_result)**2

# 結果をプロット
fig, axes = plt.subplots(nrows=2)
axes[0].plot(t, signal)
axes[0].set(title=r'Signal $f(t)$', xlabel=r"$t$", ylabel=r"f(t)")
axes[1].stem(frequencies[:len(frequencies)], power_spectrum[:len(frequencies)])
# axes[1].stem(frequencies[:len(frequencies)//2], power_spectrum[:len(frequencies)//2])
axes[1].set(title=f'Power Spectrum $P(\omega)$', xlabel=r'Frequency $ω$', ylabel=r'Power Spectrum $P(ω)$')
fig.tight_layout()
fig.show()
../../_images/5fff85839a3578e8eccedba78f57c0e98bff649c27c5e9045ad4731235777518.png

複素計量空間（ユニタリ空間）とのつながり#

信号f,gの内積を

(f,g):=∫−∞∞f(t)g(t)― dt

と定義すれば、ノルムも

‖f‖=(f,f)

と定義でき、フーリエ変換も

(F,G):=12π∫−∞∞F(ω)G(ω)―dω

とおけば

‖F‖=12π∫−∞∞|F(ω)|2 dω

となり、パーセバルの式は

(f,g)=(F,G),‖f‖2=‖F‖2

とシンプルでエレガントな表記にすることができる。

この内積は計量ベクトル空間の内積の公理のうち対称性(f,g)=(g,f)を満たさないが、ユニタリ空間におけるエルミート対称性(f,g)=(g,f)―は満たす。つまり、ユニタリ空間上で考えるなら上のシンプルな表記が可能。

（エルミート）内積

スカラとして複素数を取る線形空間Lを 複素線形空間 といい、任意の元u,v∈Lについて次の性質を満たす複素数(u,v)が定義されるとき、それを （エルミート）内積 という。

  1. 任意の u∈L に対して (u,u)≥0. 等号が成り立つのは u=0 のときのみ （ 正値性 ）

  2. 任意の u,v∈L に対して (u,v)=(v,u)― （ エルミート対称性 ）

  3. 任意の複素数 c1,c2 と任意の u,v1,v2∈L に対して (c1u1+c2u2,v)=c1(u1,v)+ c2(u2,v) （ 線形性 ）

（ユニタリ）ノルム

L を複素線形空間とする。 u∈L の （ユニタリ）ノルム を次のように定義する

‖u‖=(u,u)

ユニタリ空間（複素計量空間）

上記のような内積とノルムが定義されている複素線形空間Lを ユニタリ空間 または 複素計量空間 という。

自己相関関数#

自己相関関数

次の関数 R(τ) を信号 f(t) の 自己相関関数 と呼ぶ。

R(τ)=∫−∞∞f(t)f(t−τ)― dt

R(0) は信号 f(t) のエネルギー ∫−∞∞|f(t)|2 dt に等しく、通常は大きな値をとる。 |τ| が大きくなると R(τ) は急速に減衰する。

Hide code cell source 
import numpy as np
from scipy import signal
rng = np.random.default_rng()
sig = rng.standard_normal(100)
autocorr = signal.fftconvolve(sig, sig[::-1], mode='full')

import matplotlib.pyplot as plt
fig, (ax_orig, ax_mag) = plt.subplots(2, 1, dpi=72, figsize=[6, 4])
ax_orig.plot(sig)
ax_orig.set(title=f'White noise $f(t)$', xlabel=r"$t$")
ax_mag.plot(np.arange(-len(sig)+1,len(sig)), autocorr)
ax_mag.set(title=r'Autocorrelation $R(\tau)$', xlabel=r"$\tau$")
fig.tight_layout()
fig.show()
../../_images/61fd7aa6b7b1c6e788660a6f27c59b955fa0080de58a3d13ae69a13e8e3b3296.png
なぜ自己相関関数と呼ぶのか

実数の信号 f(t),g(t) の内積 ∫−∞∞f(t)g(t)dt の場合、f(t)とg(t)が無関係ならf(t)g(t)はそれぞれプラスとマイナスで打ち消し合ったりして小さい値になるが、f(t)とg(t)が似てるならf(t)g(t)≈f(t)2≥0となり、積分すると大きな値になる。

複素数の信号の場合はエルミート内積∫−∞∞f(t)g(t)―dtを考えると、もしf(t)とg(t)が似てるならf(t)g(t)≈f(t)f(t)―≈|f(t)|2≥0となり、積分すると大きな値になる。

このように、内積・エルミート内積は類似度の尺度となり、 （相互）相関 と呼ばれる。

f(t) と、それをτだけ平行移動した f(t−τ) の類似度なので自己相関と呼ばれる。

自己相関関数の性質
R(−τ)=R(τ)―,|R(−τ)|=|R(τ)|
証明

第1式は

R(−τ)=∫−∞∞f(t)f(t+τ)―dt=∫−∞∞f(t′−τ)f(t′)―dt′=∫−∞∞f(t′)f(t′−τ)―dt′―=R(τ)―

で示される。また共役複素数をとっても絶対値は変化しないため第2式が成り立つ。

ウィーナー・ヒンチンの定理#

ウィーナー・ヒンチンの定理

自己相関関数のフーリエ変換はパワースペクトルに等しい。

P(ω)=∫−∞∞R(τ)e−iωτ dτ

よってパワースペクトルを算出するには、2通りの方法がある

  1. 信号f(t)のフーリエ変換F(ω)を計算して絶対値の2乗を取る

  2. 自己相関関数R(τ)を計算してフーリエ変換する

信号を精度よく積分するのは難しいため、自己相関関数で算出するほうが実用上は便利らしい

Hide code cell source 
import numpy as np
from scipy import signal
rng = np.random.default_rng(seed=0)
sig = rng.standard_normal(500)

autocorr = signal.fftconvolve(sig, sig[::-1], mode='full')

import matplotlib.pyplot as plt
fig, axes = plt.subplots(3, 1, dpi=80, figsize=[6, 5])
axes[0].plot(sig)
axes[0].set(title=f'White noise $f(t)$', xlabel=r"$t$")

n = len(sig)
x = np.arange(-n+1, n)
axes[1].plot(x, autocorr)
axes[1].set(title=r'Autocorrelation $R(\tau)$', xlabel=r"$\tau$")

# Power spectrum
omega, P = signal.periodogram(sig, window='flattop', scaling='spectrum')
axes[2].semilogy(omega, P)
axes[2].set(title=r'Power spectrum $P(\omega)$', xlabel=r"$\omega$")

fig.tight_layout()
fig.show()
../../_images/428d49b490f70528aabaed2e7bf8752dd467b50e713ab06ae1abb6d0bca9e17d.png
Hide code cell source 
from scipy import signal, fft
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, dpi=80, figsize=[6, 4])
axes[0].plot(sig)
axes[0].set(title=f'White noise $f(t)$', xlabel=r"$t$")

n = len(sig)
T = 1
F = fft.fft(sig, len(sig))[:n//2]
omega = fft.fftfreq(n, T)[:n//2]
P = np.abs(F) ** 2

# Power spectrum from F(ω)
axes[1].semilogy(omega, P)
axes[1].set(title=r'Power spectrum $P(\omega)$ from $F(\omega)$', xlabel=r"$\omega$")

fig.tight_layout()
fig.show()
../../_images/547ede5624c564788755007eed71f87be8124d98b7fc3ce7014b259c1988a48f.png

それっぽいやつやってみたがうまくPが一致しない

そもそもホワイトノイズのPは1のはず…

サンプリング定理#

信号 f(t) の離散的な点 {tk},k=0,±1,±2,±3,… での値 {fk} を サンプル点 {tk} での サンプル値 と呼ぶ。 サンプル点 {tk} の間隔 τ を サンプル間隔 と呼ぶ。サンプル値 {fk} のみから連続関数を再現することを 補間 という。

次の性質を持つ関数 ϕ(t) をサンプル間隔 τ の 補間関数 と呼ぶ

ϕ(t)={1t=00t=±τ,±2τ,±3τ,…

補間関数の例はsinc（cardinal sine）関数

sinc(x)=sin⁡πxπx
Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt

def sinc(x):
    return np.sin(np.pi * x) / (np.pi * x)

fig, ax = plt.subplots(figsize=[4,2])

x = np.linspace(-3*np.pi, 3*np.pi, 100)
ax.plot(x, sinc(x))
ax.set(title="sinc(x)", xlabel="x", ylabel="sinc(x)")
fig.show()
../../_images/453dfb6ddf58540fec5dea5a0480613200e0967f42571f1d6d1c7554496d2b95.png
f^(t)=∑k=−∞∞fkϕ(t−tk)

はすべてのサンプル点tiでサンプルの値fiをとる

定義：帯域制限

信号f(t)のフーリエ変換をF(ω)とするとき、

|ω|≥W⟹F(ω)=0

であれば、信号f(t)は 帯域幅 W に 帯域制限 されているという

サンプリング間隔τを半周期（波1つ分。sinならπの幅）とする振動は各周波数W=π/τを持つ

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=[4,2])
x = np.linspace(-2*np.pi, 2*np.pi, 100)
ax.plot(x, np.sin(x))
ax.set(title="sin(x)", xlabel="x", ylabel="sin(x)")
ax.axhline(0, color="black", alpha=0.7)

ax.axvline(0, color="darkorange", linestyle="--", alpha=0.7)
ax.axvline(np.pi, color="darkorange", linestyle="--", alpha=0.7)
ax.text(np.pi / 2, -0.5, r"$\tau$", color="darkorange", ha="center")

fig.show()
../../_images/cf3b0046298433436b40b5f7a71e8da92744b99df4d0adc9cd595661b8aba6a6.png

サンプリング定理

帯域幅 W に帯域制限された信号 f(t) はサンプル間隔

τ=πW

のサンプル点 {tk} でのサンプル値 {fk} のみから次のように再現される

f(t)=∑k=−∞∞fksinct−tkτ
証明

信号f(t)のフーリエ変換は

f(t)=12π∫−∞∞F(ω)eiωt dω,F(ω)=∫−∞∞f(t)e−iωt dt

となる。F(ω)は区間[−W,W]でフーリエ級数に展開できる

F(ω)=∑k=−∞∞Ckeiπkω/W,Ck=12W∫−WWF(ω)e−iπkω/W dω

区間 [−W,W] の外では F(ω)=0 であるから、 Ck は次のように計算できる。

Ck=12W∫−∞∞F(ω)e−iπkω/W dω=πW12π∫−∞∞F(ω)eiω(−πk/W)dω=πWf(−πkW)=τf(−kτ)=τf−k

よってF(ω)は区間[−W,W]で

F(ω)=∑k=−∞∞τf−keiπkω/W=τ∑k=−∞∞fke−iπkω/W=τ∑k=−∞∞fke−ikτω=τ∑k=−∞∞fke−itkω

となる。

まとめると、

f(t)=12π∫−WWF(ω)eiωt dω=12π∫−WW(τ∑k=−∞∞fke−itkω)eiωt dω=τ2π∑k=−∞∞fk∫−WWei(t−tk)ωdω=τ2π∑k=−∞∞fk∫−WW(cos⁡(t−tk)ω+isin⁡(t−tk)ω)dω=τπ∑k=−∞∞fk∫0Wcos⁡(t−tk)ωdω=1W∑k=−∞∞fk∫0Wcos⁡(t−tk)ωdω(∵τ=π/W)=1W∑k=−∞∞fk[sin⁡(t−tk)ω(t−tk)]0W=∑k=−∞∞fksin⁡W(t−tk)W(t−tk)=∑k=−∞∞fksincW(t−tk)π=∑k=−∞∞fksinct−tkτ

一定の条件を満たす間隔で得られたサンプルと、一定の条件を満たす信号は、 離散的なサンプルからもとの連続的な信号を再現できる というのがサンプリング定理。

音声などの信号に限らず、デジタル画像などアナログ（連続値）なものをデジタル（離散値）にサンプリングする分野で広く応用されている。

間隔τでサンプリングするとW=π/τ以上の周波数の信号は捉えられない。W=π/τをサンプリング間隔τに対する ナイキスト周波数 とよぶ。

Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt

def approx(x, tau):
    # tau: サンプリング周期
    # 間隔を指定してのサンプリングが難しかったのでarangeで生成
    x_samples = np.arange(-3 * np.pi, 3 * np.pi, tau)
    y_samples = np.sin(x_samples)
    
    y_hat = []
    for i in range(len(x)):
        y_hat_i = 0
        for j in range(len(y_samples)):
            y_hat_i += y_samples[j] * np.sinc((x[i] - x_samples[j]) / tau)
        y_hat.append(y_hat_i)
    return x_samples, y_samples, np.array(y_hat)

taus = [
    ("\pi", np.pi),
    ("\pi/1.1", np.pi/1.1),
    ("\pi/2", np.pi/2),
]
fig, axes = plt.subplots(figsize=[12, 2.5], ncols=len(taus), sharey=True)
for i, ax in enumerate(axes):
    # true signal
    x = np.arange(-3 * np.pi, 3 * np.pi, 0.1)
    y = np.sin(x)
    ax.plot(x, y, label=r"true signal $f(t)$", color="black", alpha=0.7)

    # approximate
    tau_str, tau = taus[i]
    x_samples, y_samples, y_hat = approx(x, tau=tau)
    ax.plot(x_samples, y_samples, label=r"samples", alpha=0.8, linestyle=":", marker=".")
    ax.plot(x, y_hat, label=rf"approximation ($\tau = {tau_str}$)", alpha=0.8, linestyle="--")
    ax.set(xlabel=r"$t$", title=rf"Approximation ($\tau = {tau_str} = {tau:.2f}$)")#, ylim=(-1.2, 3.2))
    ax.legend()
    if i == 0:
        ax.set(ylabel=r"$f(t)$")

fig.show()
../../_images/63ca77c19fa4e7b40a435168a54fd7ed7807c458bf42765d1c14a396b9eb8c77.png

周期TがT=2τなら、周波数がω=2π/T=π/τ=W

1周期がT=2τ=2πのsinに対してサンプリング間隔がτ=πだとナイキスト周波数であり、もとの信号が復元できない様子。ω=2π/T=2π/2π=1…？

参考#

Contents