Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

フーリエ解析

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

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

フーリエ級数

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

f(t)=a02+∑k=1∞(akcos⁡kωot+bksin⁡kωot),−T2≤t≤T2f(t)=\frac{a_0}{2}+\sum_{k=1}^{\infty}\left(a_k \cos k \omega_o t+b_k \sin k \omega_o t\right), \quad-\frac{T}{2} \leq t \leq \frac{T}{2}
(1)

フーリエ係数は

ak=2T∫−T/2T/2f(t)cos⁡kωot dt,bk=2T∫−T/2T/2f(t)sin⁡kωot dta_k=\frac{2}{T} \int_{-T / 2}^{T / 2} f(t) \cos k \omega_o t \mathrm{~d} t, \quad b_k=\frac{2}{T} \int_{-T / 2}^{T / 2} f(t) \sin k \omega_o t \mathrm{~d} t
(2)

である。 ここで

ωo=2πT\omega_o=\frac{2 \pi}{T}
(3)

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

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

複素数の指数関数

複素数

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

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

xx軸を実部 x=Re⁡zx=\operatorname{Re} z、yy軸を 虚部 y=Im⁡zy=\operatorname{Im} zにとった平面を 複素平面 と呼ぶ

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()
<Figure size 400x300 with 1 Axes>

共役複素数

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

x=z+zˉ2,y=z−zˉ2ix=\frac{z+\bar{z}}{2}, \quad y=\frac{z-\bar{z}}{2 i}
(4)

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

zzˉ=(x+iy)(x−iy)=x2+y2=∣z∣2z \bar{z}=(x+i y)(x-i y)=x^2+y^2=|z|^2
(5)

が成り立つ。

オイラーの式

指数部が虚数の指数関数 eiθe^{i \theta} を次の複素数と定義する。

eiθ=cos⁡θ+isin⁡θe^{i \theta}=\cos \theta+i \sin \theta
(6)

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

定理
eiθeiϕ=ei(θ+ϕ)e^{i \theta} e^{i \phi}=e^{i(\theta+\phi)}
(7)

が成り立つ。

証明

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

eiθeiϕ=(cos⁡θ+isin⁡θ)(cos⁡ϕ+isin⁡ϕ)=(cos⁡θcos⁡ϕ−sin⁡θsin⁡ϕ)+i(cos⁡θsin⁡ϕ+sin⁡θcos⁡ϕ)ei(θ+ϕ)=cos⁡(θ+ϕ)+isin⁡(θ+ϕ)\begin{aligned} e^{i \theta} e^{i \phi} & =(\cos \theta+i \sin \theta)(\cos \phi+i \sin \phi) \\ & =(\cos \theta \cos \phi-\sin \theta \sin \phi)+i(\cos \theta \sin \phi+\sin \theta \cos \phi)\\ e^{i(\theta+\phi)} &=\cos (\theta+\phi)+i \sin (\theta+\phi)\\ \end{aligned}
(8)

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

cos⁡(θ+ϕ)=cos⁡θcos⁡ϕ−sin⁡θsin⁡ϕsin⁡(θ+ϕ)=cos⁡θsin⁡ϕ+sin⁡θcos⁡ϕ\begin{aligned} \cos (\theta+\phi) &=\cos \theta \cos \phi-\sin \theta \sin \phi\\ \sin (\theta+\phi) &=\cos \theta \sin \phi+\sin \theta \cos \phi \end{aligned}
(9)

よりわかる

定理
cos⁡θ=eiθ+e−iθ2,sin⁡θ=eiθ−e−iθ2i\cos \theta=\frac{e^{i \theta}+e^{-i \theta}}{2} , \quad \sin \theta=\frac{e^{i \theta}-e^{-i \theta}}{2 i}
(10)

が成り立つ。

証明

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

eiθ=cos⁡θ+isin⁡θ,e−iθ=cos⁡θ−isin⁡θe^{i \theta}=\cos \theta+i \sin \theta , \quad e^{-i \theta}=\cos \theta-i \sin \theta
(11)

を使うと、

eiθ+e−iθ=cos⁡θ+isin⁡θ+cos⁡θ−isin⁡θ=2cos⁡θe^{i \theta} + e^{-i \theta} = \cos \theta + i \sin \theta + \cos \theta-i \sin \theta = 2 \cos \theta
(12)

より

cos⁡θ=eiθ+e−iθ2\cos \theta=\frac{e^{i \theta}+e^{-i \theta}}{2}
(13)

が得られる。また

eiθ−e−iθ=cos⁡θ+isin⁡θ−cos⁡θ+isin⁡θ=2isin⁡θe^{i \theta} - e^{-i \theta} = \cos \theta + i \sin \theta - \cos \theta + i \sin \theta = 2i \sin \theta
(14)

より

sin⁡θ=eiθ−e−iθ2i\sin \theta=\frac{e^{i \theta}-e^{-i \theta}}{2 i}
(15)

が得られる。

フーリエ級数の複素表示

cos⁡θ=eiθ+e−iθ2,sin⁡θ=eiθ−e−iθ2i\cos \theta=\frac{e^{i \theta}+e^{-i \theta}}{2}, \quad \sin \theta=\frac{e^{i \theta}-e^{-i \theta}}{2 i}
(16)

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

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\begin{aligned} f(t) & =\frac{a_0}{2}+\sum_{k=1}^{\infty}\left(a_k \frac{e^{i k \omega_o t}+e^{-i k \omega_o t}}{2}+b_k \frac{e^{i k \omega_o t}-e^{-i k \omega_o t}}{2 i}\right) \\ & =\frac{a_0}{2}+\frac{1}{2} \sum_{k=1}^{\infty}\left(\left(a_k-i b_k\right) e^{i k \omega_o t}+\left(a_k+i b_k\right) e^{-i k \omega_o t}\right)\\ & =\sum_{k=-\infty}^{\infty} C_k e^{i k \omega_o t} \end{aligned}
(17)

ただし

Ck={(ak−ibk)/2k>0a0/2k=0(a−k+ib−k)/2k<0C_k=\left\{\begin{array}{cc} \left(a_k-i b_k\right) / 2 & k>0 \\ a_0 / 2 & k=0 \\ \left(a_{-k}+i b_{-k}\right) / 2 & k<0 \end{array}\right.
(18)

である

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

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\begin{aligned} \frac{a_k-i b_k}{2} & =\frac{2}{T} \int_{-T / 2}^{T / 2} f(t) \frac{\cos k \omega_o t-i \sin k \omega_o t}{2} \mathrm{~d} t=\frac{1}{T} \int_{-T / 2}^{T / 2} f(t) e^{-i k \omega_o t} \mathrm{~d} t \\ \frac{a_0}{2} & =\frac{1}{T} \int_{-T / 2}^{T / 2} f(t) \mathrm{d} t \\ \frac{a_{-k}+i b_{-k}}{2} & =\frac{2}{T} \int_{-T / 2}^{T / 2} f(t) \frac{\cos \left(-k \omega_o t\right)+i \sin \left(-k \omega_o t\right)}{2} \mathrm{~d} t \\ & =\frac{1}{T} \int_{-T / 2}^{T / 2} f(t)\left(\cos k \omega_o t-i \sin k \omega_o t\right) \mathrm{d} t=\frac{1}{T} \int_{-T / 2}^{T / 2} f(t) e^{-i k \omega_o t} \mathrm{~d} t \end{aligned}
(19)
f(t)=∑k=−∞∞Ckeikωot,Ck=1T∫−T/2T/2f(t)e−ikωot dtf(t)=\sum_{k=-\infty}^{\infty} C_k e^{i k \omega_o t}, \quad C_k=\frac{1}{T} \int_{-T / 2}^{T / 2} f(t) e^{-i k \omega_o t} \mathrm{~d} t
(20)

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

info - Unknown Directive
$f(t)$は実数なのに、フーリエ係数は複素数なのは問題ないのか？

→ $C_k$と$e^{i k \omega_o t}$は共役複素数なので掛けると虚部が消えて実数になる
info - Unknown Directive
位相（波の数）は角度として表される。

特に複素平面上で単位円を考えると、位相は単位円の角度（回転）として表される
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()
<Figure size 400x300 with 1 Axes>

フーリエ変換

フーリエ級数

f(t)=∑k=−∞∞Ckeikωot,Ck=1T∫−T/2T/2f(t)e−ikωot dtf(t)=\sum_{k=-\infty}^{\infty} C_k e^{i k \omega_o t}, \quad C_k=\frac{1}{T} \int_{-T / 2}^{T / 2} f(t) e^{-i k \omega_o t} \mathrm{~d} t
(21)

において

ωk=kωo,k=0,±1,±2,…\omega_k=k \omega_o, \quad k=0, \pm 1, \pm 2, \ldots
(22)

とおき、

Δω=ωk+1−ωk=ωo=2πT,Ck=F(ωk)T\Delta \omega=\omega_{k+1}-\omega_k=\omega_o=\frac{2 \pi}{T}, \quad C_k=\frac{F\left(\omega_k\right)}{T}
(23)

と書くと、

f(t)=12π∑k=−∞∞F(ωk)eiωktΔω,F(ωk)=∫−T/2T/2f(t)e−iωkt dtf(t)=\frac{1}{2 \pi} \sum_{k=-\infty}^{\infty} F\left(\omega_k\right) e^{i \omega_k t} \Delta \omega, \quad F\left(\omega_k\right)=\int_{-T / 2}^{T / 2} f(t) e^{-i \omega_k t} \mathrm{~d} t
(24)

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

F(ω)=∫−∞∞f(t)e−iωt dtF(\omega)=\int_{-\infty}^{\infty} f(t) e^{-i \omega t} \mathrm{~d} t
(25)

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

f(t)=12π∫−∞∞F(ω)eiωt dωf(t)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega t} \mathrm{~d} \omega
(26)

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

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

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()
<Figure size 640x480 with 2 Axes>

例：矩形窓

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

w(t)={1/2W−W≤t≤W0 それ以外 w(t)=\left\{\begin{array}{cc} 1 / 2 W & -W \leq t \leq W \\ 0 & \text { それ以外 } \end{array}\right.
(29)

このフーリエ変換は

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)=sinc⁡Wπω\begin{aligned} W(\omega) &=\frac{1}{2 W} \int_{-W}^W e^{-i \omega t} \mathrm{~d} t \\ &=\frac{1}{2 W} \int_{-W}^W(\cos \omega t-i \sin \omega t) \mathrm{d} t \quad (オイラーの式より e^{-i \omega t} = \cos \omega t-i \sin \omega t)\\ &=\frac{1}{2 W} \left( \int_{-W}^W \cos \omega t ~\mathrm{d} t - \int_{-W}^W i \sin \omega t ~\mathrm{d} t \right) \quad (\because 積分の線形性) \\ &=\frac{1}{2W} \int_{-W}^W \cos \omega t ~ \mathrm{d} t \quad (\because \sinは奇関数すなわち \sin(-\theta) = -\sin(\theta) であり、-Wから+Wまで積分したら0になるため消える ) \\ &=\frac{1}{W} \int_0^W \cos \omega t \mathrm{~d} t \quad (\because \int_{-W}^W \cos \omega t ~ \mathrm{d} t = 2 \int_{0}^W \cos \omega t ~ \mathrm{d} t ) \\ &=\frac{1}{W}\left[ \frac{\sin \omega t}{\omega} \right]_0^W \quad (\because \int \cos (\omega t) d t=\frac{\sin (\omega t)}{\omega} ) \\ &=\frac{1}{W} \left( \frac{\sin \omega W}{\omega} - \frac{\sin 0}{\omega} \right) \\ &=\frac{\sin W \omega}{W \omega} \quad (\because \sin 0 = 0) \\ &=\operatorname{sinc} \frac{W}{\pi} \omega \end{aligned}
(30)

となる。ただし、

sinc⁡x:=sin⁡πxπx\operatorname{sinc} x :=\frac{\sin \pi x}{\pi x}
(31)

と定義する。

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()
<Figure size 400x300 with 1 Axes>

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

w(t)=12πσe−t2/2σ2w(t)=\frac{1}{\sqrt{2 \pi} \sigma} e^{-t^2 / 2 \sigma^2}
(32)

このフーリエ変換は

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\begin{aligned} W(\omega) &= \int_{-\infty}^{\infty} w(t) e^{-i \omega t} \mathrm{~d} t \\ &= \frac{1}{\sqrt{2 \pi} \sigma} \int_{-\infty}^{\infty} e^{-t^2 / 2 \sigma^2} e^{-i \omega t} \mathrm{~d} t \\ &= \frac{1}{\sqrt{2 \pi} \sigma} \int_{-\infty}^{\infty} e^{-( \frac{t^2}{2 \sigma^2} + i \omega t)} \mathrm{~d} t \quad ( 指数法則 a^{n}a^m=a^{n+m})\\ &= \frac{1}{\sqrt{2 \pi} \sigma} \int_{-\infty}^{\infty} e^{-\left(t+i \sigma^2 \omega\right)^2 / 2 \sigma^2-\sigma^2 \omega^2 / 2} \mathrm{~d} t \quad (平方完成したり色々いじるとこうなる) \\ &= \left(\frac{1}{\sqrt{2 \pi} \sigma} \int_{-\infty}^{\infty} e^{-\left(t+i \sigma^2 \omega\right)^2 / 2 \sigma^2} \mathrm{~d} t\right) e^{-\sigma^2 \omega^2 / 2} \end{aligned}
(33)

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

W(ω)=(12πσ∫−∞∞e−z2/2σ2 dz)e−σ2ω2/2=e−σ2ω2/2=2πσ(12πσ−1e−ω2/2σ−2)\begin{aligned} W(\omega) & =\left(\frac{1}{\sqrt{2 \pi} \sigma} \int_{-\infty}^{\infty} e^{-z^2 / 2 \sigma^2} \mathrm{~d} z\right) e^{-\sigma^2 \omega^2 / 2} \\ & =e^{-\sigma^2 \omega^2 / 2}\\ &=\frac{\sqrt{2 \pi}}{\sigma}\left(\frac{1}{\sqrt{2 \pi} \sigma^{-1}} e^{-\omega^2 / 2 \sigma^{-2}}\right) \end{aligned}
(34)

となる

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()
<Figure size 400x200 with 1 Axes>
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()
<Figure size 400x240 with 1 Axes>

畳み込み積分

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

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

  • 分配法則： f(t)∗(ag(t)+bh(t))=af(t)∗g(t)+bf(t)∗h(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)=g(t) * f(t)

  • 結合法則：(f(t)∗g(t))∗h(t)=f(t)∗(g(t)∗h(t))(f(t) * g(t)) * h(t)=f(t) *(g(t) * h(t))

証明

f(t)∗g(t)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(ω)\begin{aligned} \int_{-\infty}^{\infty} f(t) * g(t) e^{-i \omega t} \mathrm{~d} t & =\int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty} f(s) g(t-s) \mathrm{d} s\right) e^{-i \omega t} \mathrm{~d} t \\ & =\int_{-\infty}^{\infty} f(s)\left(\int_{-\infty}^{\infty} g(t-s) e^{-i \omega t} \mathrm{~d} t\right) \mathrm{d} s \\ & =\int_{-\infty}^{\infty} f(s)\left(\int_{-\infty}^{\infty} g\left(t^{\prime}\right) e^{-i \omega\left(t^{\prime}+s\right)} \mathrm{d} t^{\prime}\right) \mathrm{d} s \\ & =\int_{-\infty}^{\infty} f(s)\left(\int_{-\infty}^{\infty} g\left(t^{\prime}\right) e^{-i \omega t^{\prime}} \mathrm{d} t^{\prime}\right) e^{-i \omega s} \mathrm{~d} s \\ & =\left(\int_{-\infty}^{\infty} f(s) e^{-i \omega s} \mathrm{~d} s\right)\left(\int_{-\infty}^{\infty} g\left(t^{\prime}\right) e^{-i \omega t^{\prime}} \mathrm{d} t^{\prime}\right) \\ & =F(\omega) G(\omega) \end{aligned}
(36)

（積分の順序を入れ換え、 t′=t−st^{\prime}=t-s と変数変換し、 t=t′+st=t^{\prime}+s を代入した）

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

フィルター

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

f~(t):=12W∫−WWf(t−s)ds\tilde{f}(t) := \frac{1}{2 W} \int_{-W}^W f(t-s) \mathrm{d} s
(37)

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

f~(t)=∫−∞∞w(s)f(t−s)ds=w(t)∗f(t)\tilde{f}(t)=\int_{-\infty}^{\infty} w(s) f(t-s) \mathrm{d} s=w(t) * f(t)
(38)

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

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

F~(ω)=W(ω)F(ω)\tilde{F}(\omega) = W(\omega) F(\omega)
(39)

と書くことができる。

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

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

例

fftconvolve — SciPy v1.14.1 Manual

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

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()
<Figure size 400x300 with 1 Axes>

パワースペクトル

証明
∫−∞∞f(t)g(t)‾ dt=12π∫−∞∞F(ω)G(ω)‾ dω\int_{-\infty}^{\infty} f(t) \overline{g(t)} ~\mathrm{d} t =\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) \overline{G(\omega)} ~\mathrm{d} \omega
(43)

については、

∫−∞∞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ω\begin{aligned} \int_{-\infty}^{\infty} f(t) \overline{g(t)} \mathrm{d} t & =\int_{-\infty}^{\infty} f(t)\left(\frac{1}{2 \pi} \int_{-\infty}^{\infty} \overline{G(\omega)} e^{-i \omega t} \mathrm{~d} \omega\right) \mathrm{d} t \\ & =\frac{1}{2 \pi} \int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty} f(t) e^{-i \omega t} \mathrm{~d} t\right) \overline{G(\omega)} \mathrm{d} \omega \\ & =\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) \overline{G(\omega)} \mathrm{d} \omega \end{aligned}
(44)

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

∫−∞∞∣f(t)∣2 dt=12π∫−∞∞∣F(ω)∣2 dω\int_{-\infty}^{\infty}|f(t)|^2 \mathrm{~d} t =\frac{1}{2 \pi} \int_{-\infty}^{\infty}|F(\omega)|^2 \mathrm{~d} \omega
(45)

が得られる

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

パーセバルの第2式

∫−∞∞∣f(t)∣2 dt=12π∫−∞∞∣F(ω)∣2 dω\int_{-\infty}^{\infty}|f(t)|^2 \mathrm{~d} t =\frac{1}{2 \pi} \int_{-\infty}^{\infty} |F(\omega)|^2 \mathrm{~d} \omega
(46)

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

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

P(ω)=∣F(ω)∣2P(\omega) = |F(\omega)|^2
(47)

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

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

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

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()
<Figure size 640x480 with 2 Axes>

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

信号f,gf,gの内積を

(f,g):=∫−∞∞f(t)g(t)‾ dt(f, g) := \int_{-\infty}^{\infty} f(t) \overline{g(t)} ~\mathrm{d} t
(48)

と定義すれば、ノルムも

∥f∥=(f,f)\|f\| = \sqrt{(f,f)}
(49)

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

(F,G):=12π∫−∞∞F(ω)G(ω)‾dω(F, G):=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) \overline{G(\omega)} \mathrm{d} \omega
(50)

とおけば

∥F∥=12π∫−∞∞∣F(ω)∣2 dω\|F\|=\sqrt{\frac{1}{2 \pi} \int_{-\infty}^{\infty}|F(\omega)|^2 \mathrm{~d} \omega}
(51)

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

(f,g)=(F,G),∥f∥2=∥F∥2(f, g)=(F, G), \quad\|f\|^2=\|F\|^2
(52)

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

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

自己相関関数

R(0)R(0) は信号 f(t)f(t) のエネルギー ∫−∞∞∣f(t)∣2 dt\int_{-\infty}^{\infty}|f(t)|^2 \mathrm{~d} t に等しく、通常は大きな値をとる。 ∣τ∣|\tau| が大きくなると R(τ)R(\tau) は急速に減衰する。

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()
<Figure size 432x288 with 2 Axes>
なぜ自己相関関数と呼ぶのか

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

複素数の信号の場合はエルミート内積∫−∞∞f(t)g(t)‾dt\int_{-\infty}^{\infty} f(t) \overline{g(t)} \mathrm{d} tを考えると、もしf(t)f(t)とg(t)g(t)が似てるならf(t)g(t)≈f(t)f(t)‾≈∣f(t)∣2≥0f(t)g(t) \approx f(t)\overline{f(t)} \approx |f(t)|^2 \geq 0となり、積分すると大きな値になる。

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

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

自己相関関数の性質
R(−τ)=R(τ)‾,∣R(−τ)∣=∣R(τ)∣R(-\tau)=\overline{R(\tau)} , \quad|R(-\tau)|=|R(\tau)|
(55)
証明

第1式は

R(−τ)=∫−∞∞f(t)f(t+τ)‾dt=∫−∞∞f(t′−τ)f(t′)‾dt′=∫−∞∞f(t′)f(t′−τ)‾dt′‾=R(τ)‾\begin{aligned} R(-\tau) &=\int_{-\infty}^{\infty} f(t) \overline{f(t+\tau)} \mathrm{d} t\\ &=\int_{-\infty}^{\infty} f\left(t^{\prime}-\tau\right) \overline{f\left(t^{\prime}\right)} \mathrm{d} t^{\prime} \\ &=\overline{\int_{-\infty}^{\infty} f\left(t^{\prime}\right) \overline{f\left(t^{\prime}-\tau\right)} \mathrm{d} t^{\prime}}=\overline{R(\tau)} \end{aligned}
(56)

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

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

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

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

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

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

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()
<Figure size 480x400 with 3 Axes>
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()
<Figure size 480x320 with 2 Axes>

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

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

サンプリング定理

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

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

ϕ(t)={1t=00t=±τ,±2τ,±3τ,…\phi(t)= \begin{cases}1 & t=0 \\ 0 & t= \pm \tau, \pm 2 \tau, \pm 3 \tau, \ldots\end{cases}
(58)

補間関数の例はsinc⁡\operatorname{sinc}（cardinal sine）関数

sinc⁡(x)=sin⁡πxπx\operatorname{sinc}(x) =\frac{\sin \pi x}{\pi x}
(59)
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()
<Figure size 400x200 with 1 Axes>
f^(t)=∑k=−∞∞fkϕ(t−tk)\hat{f}(t)=\sum_{k=-\infty}^{\infty} f_k \phi\left(t-t_k\right)
(60)

はすべてのサンプル点tit_iでサンプルの値fif_iをとる

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

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()
<Figure size 400x200 with 1 Axes>
証明

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

f(t)=12π∫−∞∞F(ω)eiωt dω,F(ω)=∫−∞∞f(t)e−iωt dtf(t)=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega t} \mathrm{~d} \omega, \quad F(\omega)=\int_{-\infty}^{\infty} f(t) e^{-i \omega t} \mathrm{~d} t
(64)

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

F(ω)=∑k=−∞∞Ckeiπkω/W,Ck=12W∫−WWF(ω)e−iπkω/W dωF(\omega)=\sum_{k=-\infty}^{\infty} C_k e^{i \pi k \omega / W}, \quad C_k=\frac{1}{2 W} \int_{-W}^W F(\omega) e^{-i \pi k \omega / W} \mathrm{~d} \omega
(65)

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

Ck=12W∫−∞∞F(ω)e−iπkω/W dω=πW12π∫−∞∞F(ω)eiω(−πk/W)dω=πWf(−πkW)=τf(−kτ)=τf−k\begin{aligned} C_k & =\frac{1}{2 W} \int_{-\infty}^{\infty} F(\omega) e^{-i \pi k \omega / W} \mathrm{~d} \omega\\ &=\frac{\pi}{W} \frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) e^{i \omega(-\pi k / W)} \mathrm{d} \omega \\ & =\frac{\pi}{W} f\left(-\frac{\pi k}{W}\right)=\tau f(-k \tau)=\tau f_{-k} \end{aligned}
(66)

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

F(ω)=∑k=−∞∞τf−keiπkω/W=τ∑k=−∞∞fke−iπkω/W=τ∑k=−∞∞fke−ikτω=τ∑k=−∞∞fke−itkω\begin{aligned} F(\omega) & =\sum_{k=-\infty}^{\infty} \tau f_{-k} e^{i \pi k \omega / W}=\tau \sum_{k=-\infty}^{\infty} f_k e^{-i \pi k \omega / W} \\ & =\tau \sum_{k=-\infty}^{\infty} f_k e^{-i k \tau \omega}=\tau \sum_{k=-\infty}^{\infty} f_k e^{-i t_k \omega} \end{aligned}
(67)

となる。

まとめると、

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=−∞∞fksinc⁡W(t−tk)π=∑k=−∞∞fksinc⁡t−tkτ\begin{aligned} f(t) & =\frac{1}{2 \pi} \int_{-W}^W F(\omega) e^{i \omega t} \mathrm{~d} \omega\\ & =\frac{1}{2 \pi} \int_{-W}^W\left(\tau \sum_{k=-\infty}^{\infty} f_k e^{-i t_k \omega}\right) e^{i \omega t} \mathrm{~d} \omega \\ & =\frac{\tau}{2 \pi} \sum_{k=-\infty}^{\infty} f_k \int_{-W}^W e^{i\left(t-t_k\right) \omega} \mathrm{d} \omega \\ & =\frac{\tau}{2 \pi} \sum_{k=-\infty}^{\infty} f_k \int_{-W}^W\left(\cos \left(t-t_k\right) \omega+i \sin \left(t-t_k\right) \omega\right) \mathrm{d} \omega \\ & =\frac{\tau}{\pi} \sum_{k=-\infty}^{\infty} f_k \int_0^W \cos \left(t-t_k\right) \omega \mathrm{d} \omega \\ & =\frac{1}{W} \sum_{k=-\infty}^{\infty} f_k \int_0^W \cos \left(t-t_k\right) \omega \mathrm{d} \omega \quad (\because \tau = \pi/W) \\ & =\frac{1}{W} \sum_{k=-\infty}^{\infty} f_k\left[\frac{\sin \left(t-t_k\right) \omega}{\left(t-t_k\right)}\right]_0^W\\ & =\sum_{k=-\infty}^{\infty} f_k \frac{\sin W\left(t-t_k\right)}{W\left(t-t_k\right)} \\ & =\sum_{k=-\infty}^{\infty} f_k \operatorname{sinc} \frac{W\left(t-t_k\right)}{\pi}\\ & =\sum_{k=-\infty}^{\infty} f_k \operatorname{sinc} \frac{t-t_k}{\tau} \end{aligned}
(68)

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

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

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

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()
<Figure size 1200x250 with 3 Axes>

周期TTがT=2τT=2\tauなら、周波数がω=2π/T=π/τ=W\omega=2\pi/T=\pi/\tau=W

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

参考

数学
直交関数展開：ベクトルと関数の交差点
数学
離散フーリエ変換