POINT
- sinのFFT (DFT) と DTFT，連続フーリエ変換の結果を比較する．
- numpy.fftの使い方を整理する．
numpy.fftを正弦波で試したのでメモ．
Discrete Fourier Transform (numpy.fft) — NumPy v2.0 Manual
sinの（連続）フーリエ変換
結果と見比べるために，計算しておく．\begin{aligned}
\mathcal{F}[\sin (2\pi f_{0} t)]
& = \mathcal{F} \left[\frac{e^{i2\pi f_{0} t} - e^{-i2\pi f_{0} t}}{2i}\right] \\
& = \frac{1}{2i}[\delta(f - f_{0}) - \delta(f + f_{0})] \\
& = \frac{1}{2}[e^{-i\pi/2}\delta(f - f_{0}) + e^{i\pi/2} \delta(f + f_{0})]
\end{aligned}
\mathcal{F}[\sin (2\pi f_{0} t)]
& = \mathcal{F} \left[\frac{e^{i2\pi f_{0} t} - e^{-i2\pi f_{0} t}}{2i}\right] \\
& = \frac{1}{2i}[\delta(f - f_{0}) - \delta(f + f_{0})] \\
& = \frac{1}{2}[e^{-i\pi/2}\delta(f - f_{0}) + e^{i\pi/2} \delta(f + f_{0})]
\end{aligned}
$\cos$も同様．
\begin{aligned}
\mathcal{F}[\cos (2\pi f_{0} t)]
&= \mathcal{F} \left[\frac{e^{i2\pi f_{0} t} + e^{-i2\pi f_{0} t}}{2}\right] \\
& = \frac{1}{2}[\delta(f - f_{0}) + \delta(f + f_{0})] \\
\end{aligned}
\mathcal{F}[\cos (2\pi f_{0} t)]
&= \mathcal{F} \left[\frac{e^{i2\pi f_{0} t} + e^{-i2\pi f_{0} t}}{2}\right] \\
& = \frac{1}{2}[\delta(f - f_{0}) + \delta(f + f_{0})] \\
\end{aligned}
sinのFFT (numpy.fft)
ライブラリについて
DFTは，サンプリングした$N$個のデータを$\{a_{m}\}_{m=0}^{N-1}$周期的な離散データとみなす．つまり，$a_{m+N} = a_{m}$を満たすデータ$\{a_{m}\}_{m=-\infty}^{\infty}$からのサンプリングとして解析する．信号を$\mathrm{sig}(t)$，サンプリング周期を$T_{\mathrm{s}}$とすれば，
\begin{aligned}
a_{m} = \mathrm{sig}(mT_{\mathrm{s}})
\end{aligned}
ということである．a_{m} = \mathrm{sig}(mT_{\mathrm{s}})
\end{aligned}
numpy.fftの$N$点DFTは
\begin{aligned}
&A_{k} = \sum_{m=0}^{N-1} a_{m} e^{-i2\pi km/N}
\quad (k=0,...,N-1)
\end{aligned}
で実装されている．逆DFT（IDFT）は&A_{k} = \sum_{m=0}^{N-1} a_{m} e^{-i2\pi km/N}
\quad (k=0,...,N-1)
\end{aligned}
\begin{aligned}
&a_{m} = \frac{1}{N} \sum_{k=0}^{N-1} A_{k} e^{i2\pi km/N}
\quad (k=0,...,N-1)
\end{aligned}
となる．&a_{m} = \frac{1}{N} \sum_{k=0}^{N-1} A_{k} e^{i2\pi km/N}
\quad (k=0,...,N-1)
\end{aligned}
コード
import numpy as np import matplotlib.pyplot as plt # 正弦波信号(解析対象)のパラメータ (given) f0 = 1 # 正弦波の周波数 T0 = 1 / f0 # 正弦波の基本周期 sig = lambda t: np.sin(2 * np.pi * f0 * t) # FFTのパラメータ (上記の信号をもとに決める) fs = 1 << 4 # サンプリング周波数 N = 1 << 5 # サンプル数 (2^x) Ts = 1 / fs # サンプリング周期 df = fs / N # 周波数分解能 print('サンプル数 =', N) print('周波数分解能 =', df) # 信号データ生成 tmin = 0 # サンプリング始点 tmax = tmin + N * Ts # サンプリング終点 t = np.arange(tmin, tmax, Ts) # FFT fft_sig = np.fft.fft(sig(t), N) # FFT freq = np.fft.fftfreq(t.shape[-1], Ts) # ↑に対応する周波数 # IFFT ifft_sig = np.fft.ifft(fft_sig, N) # plot data = {'Signal': (t, sig(t)), 'FFT (Re)': (freq, fft_sig.real), 'FFT (Im)': (freq, fft_sig.imag), 'amplitude spectrum': (freq, np.abs(fft_sig)), 'power spectrum': (freq, np.abs(fft_sig) ** 2), 'phase spectrum': (freq, [0 if np.abs(data) < 1e-5 else np.angle(data) for data in fft_sig])} fig = plt.figure(figsize = (10, 20)) ax1 = fig.add_subplot(5,1,1) ax2 = fig.add_subplot(5,2,3) ax3 = fig.add_subplot(5,2,4) ax4 = fig.add_subplot(5,1,3) ax5 = fig.add_subplot(5,1,4) ax6 = fig.add_subplot(5,1,5) axes = [ax1, ax2, ax3, ax4, ax5, ax6] for ax, key in zip(axes, data): ax.plot(*data[key], marker = 'o', mfc = 'None', ls = 'None', label = key) ax.set_title(key) if key == 'Signal': ax.plot(t, ifft_sig.real, marker = 'x', mfc = 'None', ls = 'None', label = 'Signal -> FFT -> IFFT') ax.legend() continue # 計算結果と比較 ax.vlines([-f0, f0], *ax.get_ylim(), colors = 'red', ls = 'solid', lw = 0.5) if key == 'phase spectrum': ax.hlines([- np.pi / 2, np.pi / 2], *ax.get_xlim(), colors = 'red', ls = 'solid', lw = 0.5) plt.show() fig.savefig('res.pdf')
出力結果
後で見るように，パラメータを（恣意的に）上手く選んでいるため，連続フーリエ変換と合致しているように見えているだけである．本来は，連続フーリエ変換ではなくDTFTと比較しなければならない．サンプル数 = 32 周波数分解能 = 0.5
sinのDTFT・DFT (手計算)
上では連続フーリエ変換とFFTの結果を比較したが，「連続フーリエ変換」は「FFTで計算しているDFT」とは異なるものである（フーリエ変換とDFTのつながり - Notes_JP）．よって，計算が上手く行っているかという確認をするためには，FFTはDTFT（これをサンプリングしたのがDFT）と比較すべきである．$\displaystyle \sin\theta = \frac{e^{i\theta} - e^{-i\theta}}{2i}$だから，
\begin{aligned}
f_{\pm}(t) = e^{\pm i2\pi f_{0} t}
\end{aligned}
のDTFT / DFTを求めれば$\sin$のDTFT / DFTも計算できる．（矩形窓を掛けた）$f_{\pm}$のDTFTはf_{\pm}(t) = e^{\pm i2\pi f_{0} t}
\end{aligned}
\begin{aligned}
& F_{\pm}(f) \\
&= \sum_{k=0}^{N-1} f_{\pm}(k T_{\mathrm{s}}) e^{-i2\pi f k T_{\mathrm{s}} } \\
&= \sum_{k=0}^{N-1} e^{-i2\pi (f \mp f_{0}) k T_{\mathrm{s}}} \\
&=
\begin{cases}
\, N & \!\!(\theta_{\pm} / 2\pi \in \mathbb{Z}) \\
\, \displaystyle
\frac{\sin(N\theta_{\pm} / 2)}{\sin(\theta_{\pm} / 2)} e^{-i(N-1)\theta_{\pm} / 2} & \!\!(\text{otherwise})
\end{cases} \\
&\bigl(\theta_{\pm} = 2\pi (f \mp f_{0}) T_{\mathrm{s}} \bigr)
\end{aligned}
となる．この$f = nf_{\mathrm{s}} / N$における値（$n=0,1,...,N-1$）が$f_{\pm}$のDTFである（上の計算では& F_{\pm}(f) \\
&= \sum_{k=0}^{N-1} f_{\pm}(k T_{\mathrm{s}}) e^{-i2\pi f k T_{\mathrm{s}} } \\
&= \sum_{k=0}^{N-1} e^{-i2\pi (f \mp f_{0}) k T_{\mathrm{s}}} \\
&=
\begin{cases}
\, N & \!\!(\theta_{\pm} / 2\pi \in \mathbb{Z}) \\
\, \displaystyle
\frac{\sin(N\theta_{\pm} / 2)}{\sin(\theta_{\pm} / 2)} e^{-i(N-1)\theta_{\pm} / 2} & \!\!(\text{otherwise})
\end{cases} \\
&\bigl(\theta_{\pm} = 2\pi (f \mp f_{0}) T_{\mathrm{s}} \bigr)
\end{aligned}
\begin{aligned}
&\!\!\!\!\!\!\! \underline{\theta = 2 m \pi \:(m \in \mathbb{Z})} \\
&\sum_{k=0}^{N-1} e^{-ik\theta} = N \\
&\!\!\!\!\!\!\! \underline{\text{otherwise}} \\
&\sum_{k=0}^{N-1} e^{-ik\theta}
=\frac{1 - e^{-iN\theta}}{1 - e^{-i\theta}} \\
&=\frac{e^{-iN\theta / 2} (e^{iN\theta / 2} - e^{-iN\theta / 2})}{e^{-i\theta/2} (e^{i\theta/2} - e^{-i\theta/2})} \\
&=\frac{\sin(N\theta / 2)}{\sin(\theta / 2)} e^{-i(N-1)\theta / 2}
\end{aligned}
を使った）．&\!\!\!\!\!\!\! \underline{\theta = 2 m \pi \:(m \in \mathbb{Z})} \\
&\sum_{k=0}^{N-1} e^{-ik\theta} = N \\
&\!\!\!\!\!\!\! \underline{\text{otherwise}} \\
&\sum_{k=0}^{N-1} e^{-ik\theta}
=\frac{1 - e^{-iN\theta}}{1 - e^{-i\theta}} \\
&=\frac{e^{-iN\theta / 2} (e^{iN\theta / 2} - e^{-iN\theta / 2})}{e^{-i\theta/2} (e^{i\theta/2} - e^{-i\theta/2})} \\
&=\frac{\sin(N\theta / 2)}{\sin(\theta / 2)} e^{-i(N-1)\theta / 2}
\end{aligned}
よって，
\begin{aligned}
g(t)
&= \sin(2\pi f t + \phi) \\
&= \frac{e^{i(2\pi f t + \phi)} - e^{-i(2\pi f t + \phi)}}{2i}
\end{aligned}
のDFTはg(t)
&= \sin(2\pi f t + \phi) \\
&= \frac{e^{i(2\pi f t + \phi)} - e^{-i(2\pi f t + \phi)}}{2i}
\end{aligned}
\begin{aligned}
&G(n f_{\mathrm{s}}/N) \\
&=\frac{1}{2} \biggl[
F_{+}(n f_{\mathrm{s}}/N) e^{-i(\pi/2 + \phi)}
+ F_{-}(n f_{\mathrm{s}}/N) e^{i(\pi/2 - \phi)}
\biggr]
\end{aligned}
となる（$n=0,1,...,N-1$）．&G(n f_{\mathrm{s}}/N) \\
&=\frac{1}{2} \biggl[
F_{+}(n f_{\mathrm{s}}/N) e^{-i(\pi/2 + \phi)}
+ F_{-}(n f_{\mathrm{s}}/N) e^{i(\pi/2 - \phi)}
\biggr]
\end{aligned}
以下では$\phi = 0$の場合で
- FFTの振幅スペクトル＝$|G(n f_{\mathrm{s}}/N)|$
- FFTの位相スペクトル＝$\arg [G(n f_{\mathrm{s}}/N)]$
FFT(numpy.fft)とDTFT(手計算)の比較
numpy.fftで計算したFFT(DFT)と，上で手計算したDTFTを一緒にプロットする．DTFTをサンプリングしたものがDFTだから，サンプリング周波数を変えたときにDFTがどのように変わり得るかを確認できる．比較1
上で計算した$\sin$に矩形窓を掛けた場合のDTFTを計算し，上図のプロットと重ねる．まず，手計算のDTFTだけをプロットしたものを見てみる．def Fp(N, fs, f0, f): if abs((f - f0) / fs - (f - f0) // fs) < 1e-10: return N + 0j else: theta = 2 * np.pi * (f - f0) / fs psi = complex(0, - (N - 1) * theta / 2) return np.exp(psi) * np.sin(N * theta / 2) / np.sin(theta / 2) def Fm(N, fs, f0, f): if abs((f + f0) / fs - (f + f0) // fs) < 1e-10: return N + 0j else: theta = 2 * np.pi * (f + f0) / fs psi = complex(0, - (N - 1) * theta / 2) return np.exp(psi) * np.sin(N * theta / 2) / np.sin(theta / 2) def DTFT_of_sin(N, fs, f0, f): return (Fp(N, fs, f0, f) * (-1j) + Fm(N, fs, f0, f) * (1j)) / 2 Fp = np.vectorize(Fp) Fm = np.vectorize(Fm) DFT_of_sin = np.vectorize(DTFT_of_sin) # --- test --- f0 = 1 fs = 1 << 4 N = 1 << 5 freq = np.linspace(0, fs // 2, 1 << 8) DTFT = DTFT_of_sin(N, fs, f0, freq) fig = plt.figure(figsize = (10, 5)) ax1 = fig.add_subplot(2,1,1) ax2 = fig.add_subplot(2,1,2) ax1.set_title('amplitude spectrum') ax2.set_title('phase spectrum') ax1.plot(freq, np.abs(DTFT), marker = 'o', mec = 'tab:blue', mfc = 'None') ax1.plot(freq, np.abs(DTFT), marker = 'None', c = 'tab:orange') ax2.plot(freq, np.angle(DTFT), marker = 'o', mec = 'tab:blue', mfc = 'None') ax2.plot(freq, np.angle(DTFT), marker = 'None', c = 'tab:orange') plt.tight_layout() plt.show() fig.savefig('res.pdf')
これをnumpy.fftの結果と重ねてみる．すると，上の計算結果は，DTFTのうち連続フーリエ変換と整合する点のみをサンプリングしたものになっていることがわかる．
DTFTの$\theta_{\pm} / 2\pi = (f \mp f_{0}) / f_{\mathrm{s}}$は，DFTでサンプリングする点$f = nf_{\mathrm{s}} / N$（$n=0,1,...,N-1$）においては$\theta_{\pm} / 2\pi = n/N \mp f_{0} / f_{\mathrm{s}}$である．$2 f_{0} \leq f_{\mathrm{s}}$となるようにサンプリングする（ナイキスト周波数）から，$n / N = \pm f_{0} / f_{\mathrm{s}}$を満たす$n = n_{\pm}$でのみ$F_{\pm} (n_{\pm} f_{\mathrm{s}} / N) = N$が成り立つ．当然，$f_{0}, f_{\mathrm{s}}, N$のとり方によっては，この条件を満たす$n = n_{\pm}$は存在しないこともある．
上の例では，$f_{0} = 1, f_{\mathrm{s}} = 2^{4} = 16, N = 2^{5} = 32$であるから，$n_{\pm} = \pm 2$となる．FFTの振幅スペクトルのピークは$n_{\pm} = \pm 2$に現れる．
振幅スペクトルがゼロとなるのは，$F_{\pm}(f) = 0$となる場合である．これは，$\sin(N\theta_{\pm}/2) = 0$（かつ$\sin(\theta_{\pm}/2) \neq 0$）の場合に成り立つ．$f = nf_{\mathrm{s}} / N$（$n=0,1,...,N-1$）において$N\theta_{\pm}/2\pi = n \mp N f_{0} / f_{\mathrm{s}}$であるから，「正弦波の基本周波数$f_{0}$」が「周波数分解能$f_{\mathrm{s}}/N$」の倍数なら，$n \neq n_{\pm}$の点では常にゼロとなる．上の計算例では確かにそうなっている．
以上で，上の計算例は，パラメータを都合よく取ることで，DTFTから連続フーリエ変換と整合する（ように見えている）ことがわかった．
比較2：スペクトル漏れ
次に，パラメータを変えて，連続フーリエ変換と「一致しない」DTFTの曲線上にサンプリング点が乗る場合を見てみる．上の例では「周波数分解能$f_{\mathrm{s}}/N = 0.5$」であるから，例えば他の条件は変えずに「正弦波の基本周波数$f_{0} = 1.25$」とすれば$n \neq n_{\pm}$の点でゼロでない値を取る（スペクトル漏れ）．
また，図からもわかるように，そもそも$n / N = \pm f_{0} / f_{\mathrm{s}}$を満たす$n = n_{\pm}$は存在しない．
比較3：ゼロパディング
上の例では，スペクトルのピークをFFTで正しく捉えることができていない．これは，周波数分解能が$f_{\mathrm{s}} / N = 0.5$なのに対し，正弦波の基本周波数が$f_{0} = 1.25$だからである．原因は，DFTで周波数を$f_{\mathrm{s}} / N$を単位としてDTFTをサンプリングするという，DFTの枠組み（定義の仕方）による制限である．つまり，サンプリングの仕方を変えれば分解能を上げることはできる．
とはいっても，FFTはDFTの枠組みに従っているので，DFTの枠組みで分解能を向上する方法がほしい．少し考えてみれば，「窓関数を掛けたあとの信号データ」のデータ数を増やせば，周波数分解能を向上させた結果が得られることに気づく（DTFTは同じものだが，サンプリングが細かくなる）．これは，窓の外のデータを増やすことに相当し，データとしてゼロを追加して$N$を大きくすれば良いだけである．
上の例では$N = 2^{5} = 32$で$f_{\mathrm{s}} / N = 0.5$だから，$N = 2^{6} = 64$とすれば$f_{\mathrm{s}} / N = 0.25$となり，スペクトルのピークを正しく捉えられるはずである．
numpy.fft.fftにおいて，$n$としてデータ点数より大きい値を指定すれば，自動的にゼロでパディングされる（numpy.fft.fft — NumPy v2.0 Manual）．したがって，最初のコードを少し書き換えるだけで良い．
import numpy as np import matplotlib.pyplot as plt # 正弦波信号(解析対象)のパラメータ (given) f0 = 1.25 # 正弦波の周波数 T0 = 1 / f0 # 正弦波の基本周期 sig = lambda t: np.sin(2 * np.pi * f0 * t) # FFTのパラメータ (上記の信号をもとに決める) fs = 1 << 4 # サンプリング周波数 Ns = 1 << 5 # サンプル数 N = 1 << 6 # ゼロパディング後のサンプル数 (2^x, N >= Ns) Ts = 1 / fs # サンプリング周期 df = fs / N # 周波数分解能 print('サンプル数 =', Ns) print('ゼロパディング後のサンプル数 =', N) print('周波数分解能 =', df) # 信号データ生成 tmin = 0 # サンプリング始点 tmax = tmin + Ns * Ts # サンプリング終点 t = np.arange(tmin, tmax, Ts) # [start, stop) t_zero_pad = np.arange(tmin, tmin + N * Ts, Ts) # FFT fft_sig = np.fft.fft(sig(t), N) # FFT freq = np.fft.fftfreq(N, Ts) # ↑に対応する周波数 # IFFT ifft_sig = np.fft.ifft(fft_sig, N) # DTFT def Fp(Ns, fs, f0, f): if abs((f - f0) / fs - (f - f0) // fs) < 1e-10: return Ns + 0j else: theta = 2 * np.pi * (f - f0) / fs psi = complex(0, - (Ns - 1) * theta / 2) return np.exp(psi) * np.sin(Ns * theta / 2) / np.sin(theta / 2) def Fm(Ns, fs, f0, f): if abs((f + f0) / fs - (f + f0) // fs) < 1e-10: return Ns + 0j else: theta = 2 * np.pi * (f + f0) / fs psi = complex(0, - (Ns - 1) * theta / 2) return np.exp(psi) * np.sin(Ns * theta / 2) / np.sin(theta / 2) def DTFT_of_sin(Ns, fs, f0, f): return (Fp(Ns, fs, f0, f) * (-1j) + Fm(Ns, fs, f0, f) * (1j)) / 2 Fp = np.vectorize(Fp) Fm = np.vectorize(Fm) DFT_of_sin = np.vectorize(DTFT_of_sin) DTFT_freq = np.linspace(0, fs // 2, 1 << 10) DTFT = DFT_of_sin(Ns, fs, f0, DTFT_freq) DTFT_amp = np.abs(DTFT) DTFT_phase = np.angle(DTFT) # plot data = {'Signal': (t, sig(t)), 'FFT (Re)': (freq, fft_sig.real), 'FFT (Im)': (freq, fft_sig.imag), 'amplitude spectrum': (freq, np.abs(fft_sig)), 'power spectrum': (freq, np.abs(fft_sig) ** 2), 'phase spectrum': (freq, [0 if np.abs(data) < 1e-5 else np.angle(data) for data in fft_sig])} fig = plt.figure(figsize = (10, 20)) ax1 = fig.add_subplot(5,1,1) ax2 = fig.add_subplot(5,2,3) ax3 = fig.add_subplot(5,2,4) ax4 = fig.add_subplot(5,1,3) ax5 = fig.add_subplot(5,1,4) ax6 = fig.add_subplot(5,1,5) axes = [ax1, ax2, ax3, ax4, ax5, ax6] for ax, key in zip(axes, data): ax.plot(*data[key], marker = 'o', mfc = 'None', ls = 'None', label = key, zorder = 10) ax.set_title(key) if key == 'Signal': ax.plot(t_zero_pad, ifft_sig.real, marker = 'x', mfc = 'None', ls = 'None', label = 'Signal -> FFT -> IFFT') ax.legend() continue # 計算結果と比較 if key == 'amplitude spectrum': ax.plot(DTFT_freq, DTFT_amp, zorder = 1) ax.vlines([-f0, f0], *ax.get_ylim(), colors = 'red', ls = 'solid', lw = 0.5) if key == 'phase spectrum': ax.hlines([- np.pi / 2, np.pi / 2], *ax.get_xlim(), colors = 'red', ls = 'solid', lw = 0.5) ax.plot(DTFT_freq, DTFT_phase, zorder = 1) plt.show() fig.savefig('res.pdf')
サンプル数 = 32 ゼロパディング後のサンプル数 = 64 周波数分解能 = 0.25
参考文献
- [1] 信号とシステム (斉藤洋一)：7.3. 信号のスペクトル解析
- [2] Pythonで学ぶフーリエ解析と信号処理 (神永正博)
