離散フーリエ変換

離散フーリエ変換#

離散フーリエ変換#

複素フーリエ級数#

θを角度とすると、円周に沿って値が定義された関数f(θ)（sinやcosなどのこと？？）は周期T=2πの周期関数であり、θに2πの任意の整数倍を足しても引いてもf(θ)は同じ値になる。

フーリエ係数の複素表示

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

より、f(θ)の基本周波数はω0=2π/T=1であるから、f(θ)のフーリエ級数は、

f(θ)=∑k=−∞∞Ckeikθ,Ck=12π∫−ππf(θ)e−ikθdθ

と書くことができる。

離散フーリエ変換#

円周上をN分割し、N個のサンプル点をとる

θl=2πNl,l=0,1,2,…,N−1

（1周期が2πなのをN分割したもののl倍がθl）

このサンプル点での f(θ) のサンプル値を fl=f(θl) とする。

前述のフーリエ係数は連続関数 f(θ) を無限個の係数 {Ck},k=0,±1,±2,±3,…, で表すものだが、 もし N 個のサンプル値 {fl} のみが必要な場合は N 個の係数のみで表される。

fl=∑k=0N−1Fkei2πkl/N,Fk=1N∑l=0N−1fle−i2πkl/N

係数 {Fk} をデータ {fl} の 離散フーリエ変換 と呼ぶ。

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_l = 2\pi l / N$")

# 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/6901df82e8988383e1fd6667cf90ee27d3778eb4dd6ccae95c4777a14df59302.png

Tip

書籍によっては1/Nの表記方法が異なる場合がある。1/Nをflのほうにもってきて

fl=1N∑k=0N−1Fkei2πkl/N,Fk=∑l=0N−1fle−i2πkl/N

としたり、あるいはl,kについて平等にするため

fl=1N∑k=0N−1Fkei2πkl/N,Fk=1N∑l=0N−1fle−i2πkl/N

とすることがある

逆フーリエ変換#

クロネッカーのデルタ関数の離散バージョン
1N∑k=0N−1ei2π(m−n)k/N={1m≡n(modN)0m≢n(modN)

ただし m≡n(modN) （N を 法 として 合同 であると読む）は m−n が N の倍数であることを表す。

これを使うことで、データ{fl}から離散フーリエ変換Fkを定義すると

∑k=0N−1Fkei2πkl/N=∑k=0N−1(1N∑m=0N−1fme−i2πkm/N)ei2πkl/N=∑m=0N−1fm(1N∑k=0N−1ei2π(l−m)k/N)

となる。最後の項のカッコの中はl≡m (modN)のとき1、それ以外は0となる。0≤l<N,0≤m<N の範囲では l≡m(modN) となるのはl=mの場合のみなので、fmを掛けて和∑m=0N−1をとるとflになる。よって逆フーリエ変換の式fl=∑k=0N−1Fkei2πkl/Nが成立する。

周期的な添字に拡張する#

取り扱いを便利にするため、以下ではfl,Fkのl,k=0,1,…,N−1の値を周期的に拡張する。 例えばfN=f0,fN+1=f1,…とする。

このように拡張すると、総和は任意の連続する N 個の和に置き換えても同じになる。例えば ∑k=0N−1 は ∑k=1N,∑k=2N+1,∑k=3N+2,… と書いても ∑k=−1N−2,∑k=−2N−3,… と書いても同じである。

周期関数のサンプリング定理#

帯域制限#

周期 2π の連続関数 f(θ) がフーリエ級数に展開されるとき、そのフーリエ係数 Ck がある k の範囲以外は 0 であるなら f(θ) は 帯域制限 されているという。

帯域制限された周期関数は、ある間隔より細かくサンプルすればフーリエ係数 Ck と離散フーリエ変換 Fk が等しくなる。

離散フーリエ変換Fk,|k|<N2は次のように書くことができる。

Fk=1N∑l=0N−1f(2πlN)e−i2πkl/N=1N∑l=0N−1(∑m=−∞∞Cmei2πlm/N)e−i2πkl/N=∑−N/2<m<N/2Cm(1N∑l=0N−1ei2π(m−k)l/N)

−N/2<m<N/2,−N/2<k<N/2 のとき m≡k (modN) となるのは m=k の場合しかない。 ゆえに 上式は Ck に等しい。

よって、 帯域制限された周期関数は、ある間隔より細かくサンプルすればそのサンプル値の補間によって表現できる。

周期関数のサンプリング定理#

周期関数のサンプリング定理

周期 2π の連続関数 f(θ) のフーリエ係数 Ck が |k|≥N2 に対して 0 のとき、f(θ) は区間 [0,2π] を N 等分して得られるサンプル値 fl=f(θl) から次のように再現される

f(θ)=∑l=0N−1flϕN(θ−θl)

ただし、 ϕN(θ) は次のように定義した補間関数である。

ϕN(θ)=1N∑−N/2<k<N/2eikθ=1+2∑0<k<N/2cos⁡kθN

|k|≥N/2 では Ck=0 であり、 |k|< N/2 では Ck=Fk であるから、 f(θ) は次のように書ける。

f(θ)=∑−N/2<k<N/2Fkeikθ=∑−N/2<k<N/2(1N∑l=0N−1fle−i2πkl/N)eikθ=∑l=0N−1fl(1N∑−N/2<k<N/2eik(θ−2πl/N))=∑l=0N−1fl(1N∑−N/2<k<N/2eik(θ−θl))=∑l=0N−1flϕN(θ−θl)

畳み込み和定理#

定義：畳み込み和

周期 N のデータ {fl},{gl} の（循環）たたみこみ和 {fl∗gl} を次のように定義する。

fl∗gl=1N∑m=0N−1fmgl−m

畳み込み和定理

周期 N のデータ {fl},{gl} の離散フーリエ変換をそれぞれ {Fk},{Gk} とするとき、 {fl∗gl} の離散フーリエ変換は {FkGk} である

証明
1N∑l=0N−1fl∗gle−i2πkl/N=1N∑l=0N−1(1N∑m=0N−1fmgl−m)e−i2πkl/N=1N∑m=0N−1fm(1N∑l=0N−1gl−me−i2πkl/N)=1N∑m=0N−1fm(1N∑l′=−mN−1−mgl′e−i2πk(l′+m)/N)=1N∑m=0N−1fm(1N∑l′=0N−1gl′e−i2πkl′/N)e−i2πkm/N=(1N∑m=0N−1fme−i2πkm/N)(1N∑l′=0N−1gl′e−i2πkl′/N)=FkGk

パワースペクトル#

パーセバルの式

1N∑l=0N−1flgl―=∑k=0N−1FkGk―,1N∑l=0N−1|fl|2=∑k=0N−1|Fk|2

(1/N)∑l=0N−1|fl|2 は {fl} の平均エネルギーを表しているとみなせる。パーセバルの式の第2式はこれが ∑k=0N−1|Fk|2 で表されることを意味している。したがってパワースペクトルを

Pk:=|Fk|2

と定義するとパーセバルの式の第2式は

1N∑l=0N−1|fl|2=∑k=0N−1Pk

と書き換えることができる。

データ{fk}が実数のとき、F−k=Fk―であるからP−k=Pkである。なのでPkのグラフをk=…,−2,−1,0,1,2,…に対してプロットするとk=0に関して左右対称になる。

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) # 時間

f = np.zeros_like(t)
components = [
    # (スペクトルF_k, 周期θ)
    (1.5, 50),
    (2, 100),
    (2.5, 200),
    (2, 300),
    (1, 500),
]
for F_k, frequency_k in components:
    f += F_k * np.sin(t * 2 * np.pi * frequency_k)

F = np.fft.fft(f) # フーリエ変換
freq = np.fft.fftfreq(N, d=d)  # 周波数スケール

# 振幅スペクトル Amplitude の取得
# 振幅スペクトルは信号をフーリエ変換した結果の絶対値をとったもの
Amp = np.abs(F)
Amp = Amp / (N / 2) # 正規化
Pow = Amp**2

# 結果をプロット
fig, axes = plt.subplots(nrows=2)
axes[0].plot(t, f)
axes[0].set(title=r'Data $\{ f_l \}$', xlabel=r"$l$", ylabel=r"$f_l$")

axes[1].stem(freq, Pow, 'r', markerfmt=" ", basefmt=" ")
# 左右対称なので正の値だけ N//2 で取り出してもいい
# axes[1].stem(freq[:N//2], Amp[:N//2], 'r', markerfmt=" ", basefmt=" ")
axes[1].set(title=f'Power $P_k$', xlabel=r'$k$', ylabel=r'$P_k$', xlim=(-550, 550))

fig.tight_layout()
fig.show()
../../_images/55d25c806f3dbaac37d056cfcb9d4b542687420c35e07c02204d4e6c5a7d3e2e.png
Hide code cell source 
import numpy as np
import matplotlib.pyplot as plt

# サンプルデータを生成する
N = 100
d = 0.001 # サンプリング周期
t = np.arange(0, N * d, d) # 時間

f = np.zeros_like(t)
for k in range(2 * N):  # 適当にNより大きいkにしてみる
    F_k = k
    frequency_k = k
    f += F_k * np.sin(t * 2 * np.pi * frequency_k)

F = np.fft.fft(f) # フーリエ変換
freq = np.fft.fftfreq(N, d=d)  # 周波数スケール

# 振幅スペクトル Amplitude の取得
# 振幅スペクトルは信号をフーリエ変換した結果の絶対値をとったもの
Amp = np.abs(F)
Amp = Amp / (N / 2) # 正規化
Pow = Amp**2

# 結果をプロット
fig, axes = plt.subplots(nrows=2)
axes[0].plot(t, f)
axes[0].set(title=r'Data $\{ f_l \}$', xlabel=r"$l$", ylabel=r"$f_l$")

axes[1].stem(freq, Pow, 'r', markerfmt=" ", basefmt=" ")
# 左右対称なので正の値だけ N//2 で取り出してもいい
# axes[1].stem(freq[:N//2], Amp[:N//2], 'r', markerfmt=" ", basefmt=" ")
axes[1].set(title=f'Power $P_k$', xlabel=r'$k$', ylabel=r'$P_k$')

fig.tight_layout()
fig.show()
../../_images/ecaf55413f17d45771ace776ffce79a958e89b135df73a95d5622e8548068137.png

自己相関係数#

周期 N のデータ {fl} の自己相関係数 {Rn} を次のように定義する。

Rn=1N∑l=0N−1flfl−n―

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

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

周期 N のデータ {fl} の自己相関係数 {Rn} の離散フーリエ変換はパワースペクトルに等しい

Pk=1N∑n=0N−1Rne−i2πkn/N
証明
1N∑n=0N−1Rne−i2πkn/N=1N∑n=0N−1(1N∑l=0N−1flfl−n―)e−i2πkn/N=1N∑l=0N−1fl(1N∑n=0N−1fl−n―e−i2πkn/N)=1N∑l=0N−1fl(1N∑n=0N−1fl−nei2πkn/N)―=1N∑l=0N−1fl(1N∑n′=0N−1fn′ei2πk(l−n′)/N)=(1N∑l=0N−1fle−i2πkl/N)(1N∑n′=0N−1fn′e−i2πkn′/N)―=FkFk―=|1N∑n′=0N−1f2=Pk

まとめ#

fl⟶Fk=1N∑l=0N−1fle−i2πkl/N↓↓Rn=1N∑l=0N−1flfl−n―⟶Pk={|Fk|21N∑n=0N−1Rne−i2πkn/N

連続のフーリエ変換と同様に、ウィーナー・ヒンチンの定理を使うことでもパワースペクトルを得ることができる。

しかし離散の場合は高速フーリエ変換があるため、ウィーナー・ヒンチンの定理を使うことによる計算量削減などの効果は相対的に低い。実用上は高速フーリエ変換一択になる。

Contents