NumPy 逆フーリエ変換(IDFT):DFTからの信号再構成プログラミング
NumPyにおいて、離散フーリエ変換(DFT)は、時間や空間領域で表現された離散的なデータを、周波数領域の成分に分解するための重要な数学的な操作です。NumPyのnumpy.fft
モジュールに実装されており、信号処理、画像処理、データ解析など、さまざまな分野で活用されています。
DFTの基本的な考え方
私たちが普段扱う多くのデータは、時間の経過に伴う変化(例えば、音の波形、株価の変動など)や、空間的な分布(例えば、画像のピクセルの明るさ)として表現されます。DFTは、これらのデータを「どのような周波数の波が、どれくらいの強さで含まれているか」という情報に変換します。
例えるなら、音楽を聴いたときに、その音楽が様々な高さの音(周波数)の組み合わせでできていると感じるのと同じです。DFTは、その音の波形データから、それぞれの周波数の成分を数値的に抽出するようなイメージです。
NumPyにおけるDFTの関数
NumPyのnumpy.fft
モジュールでよく使われるDFT関連の関数は主に以下のものです。
-
numpy.fft.ifftshift(x, axes=None)
fftshift()
で並べ替えられた配列を元の順序に戻します。 -
numpy.fft.fftshift(x, axes=None)
DFTの結果の直流成分(0 Hz)が配列の先頭に来るように並べ替えられた配列を、中央に移動します。これにより、周波数スペクトルを視覚的に理解しやすくなります。 -
numpy.fft.fftfreq(n, d=1.0)
DFTの結果として得られる各周波数成分に対応する周波数の値の配列を生成します。n
: DFTの点の数。d
: サンプリング間隔(時間の単位)。例えば、1秒間に100回サンプリングした場合、d = 1/100
となります。
-
numpy.fft.ifft(a, n=None, axis=-1, norm=None)
numpy.fft.fft()
で得られた周波数領域のデータa
に対して、逆離散フーリエ変換を実行し、時間や空間領域のデータに戻します。引数はfft()
と同様です。 -
numpy.fft.fft(a, n=None, axis=-1, norm=None)
1次元の複素数列a
に対して、離散フーリエ変換を実行します。a
: 変換したい入力配列(通常は実数または複素数のデータ)。n
: 変換する点の数(入力配列の長さより短い場合は切り捨て、長い場合はゼロパディングされます)。省略した場合は入力配列の長さが使われます。axis
: 変換を行う軸を指定します。省略した場合は最後の軸が使われます。norm
: 正規化の方法を指定します。"ortho"
を指定すると、順変換と逆変換でユニタリ性(エネルギー保存)が保たれるように正規化されます。デフォルトはNone
です。
NumPyでDFTを使う基本的な流れ
- データの準備
時間領域や空間領域の離散的なデータ(NumPy配列)を用意します。 - DFTの実行
numpy.fft.fft()
関数を使って、データを周波数領域のデータに変換します。 - 周波数情報の解析
得られた周波数成分の大きさ(絶対値)や位相などを解析します。numpy.abs()
で振幅スペクトル、numpy.angle()
で位相スペクトルを得ることができます。 - (必要に応じて) 逆DFTの実行
numpy.fft.ifft()
関数を使って、周波数領域のデータを元の時間領域や空間領域のデータに戻します。
簡単な例
import numpy as np
import matplotlib.pyplot as plt
# サンプリング周波数と時間
fs = 100 # サンプル/秒
t = np.linspace(0, 1, fs, endpoint=False)
# 2Hzと10Hzのサイン波を合成した信号
f1 = 2
f2 = 10
signal = np.sin(2 * np.pi * f1 * t) + 0.5 * np.sin(2 * np.pi * f2 * t)
# DFTを実行
fft_result = np.fft.fft(signal)
# 周波数軸の作成
freq = np.fft.fftfreq(signal.size, d=1/fs)
# 振幅スペクトルを計算
amplitude_spectrum = np.abs(fft_result)
# プロット
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title("時間領域の信号")
plt.xlabel("時間 [秒]")
plt.ylabel("振幅")
plt.subplot(2, 1, 2)
plt.plot(freq[:signal.size//2], amplitude_spectrum[:signal.size//2]) # ナイキスト周波数までを表示
plt.title("周波数領域の信号 (振幅スペクトル)")
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅")
plt.grid(True)
plt.tight_layout()
plt.show()
配列の形状(Shape)に関するエラー
- トラブルシューティング
- 入力配列の形状を確認し、意図した次元数であることを確かめます。
.shape
属性で確認できます。 - 多次元配列に対してDFTを行いたい場合は、
axis
引数を適切に設定します。例えば、2次元の画像データに対して行ごとにDFTを行う場合はaxis=1
、列ごとに行う場合はaxis=0
を指定します。
- 入力配列の形状を確認し、意図した次元数であることを確かめます。
- 原因
numpy.fft.fft()
などの関数は、通常1次元の配列に対して動作します。多次元配列を渡す場合は、どの軸に対してDFTを行うかをaxis
引数で明示する必要があります。 - エラー
ValueError: object too deep for desired array
や、DFT関数の引数として不適切な形状の配列を渡した場合に発生することがあります。
変換長の誤り(Incorrect Transform Length)
- トラブルシューティング
- 通常は、入力データの長さと同じ
n
の値を使用することが推奨されます。 - 特定の周波数分解能が必要な場合や、異なる長さのデータ間で比較を行いたい場合にのみ、
n
の値を調整することを検討してください。その際も、ゼロパディングが結果に与える影響を理解しておく必要があります。
- 通常は、入力データの長さと同じ
- 原因
n
が入力データの長さより短い場合、データが切り捨てられて周波数分解能が低下します。n
が入力データの長さより長い場合、データにゼロがパディングされ、周波数成分の間に不要な点が挿入される可能性があります。
- 問題
n
引数を使用してDFTの変換長を指定する際に、入力データの長さと異なる値を指定すると、意図しない結果になることがあります。
正規化(Normalization)に関する誤解
- トラブルシューティング
- 順変換後の結果を逆変換して元の信号に戻したい場合は、順変換と逆変換で同じ
norm
の値を指定するか、適切にスケールを調整する必要があります。 - エネルギー保存を重視する場合は
norm="ortho"
を使用することを検討してください。この場合、順変換と逆変換の両方で N​1(Nは変換長) で正規化されます。
- 順変換後の結果を逆変換して元の信号に戻したい場合は、順変換と逆変換で同じ
- 原因
デフォルトではnorm=None
であり、この場合、DFTとIDFTは互いにスケールが異なります。エネルギー保存則を満たすような正規化を行いたい場合はnorm="ortho"
を指定する必要があります。 - 問題
norm
引数の使い方を誤ると、順変換と逆変換の結果のスケールが一致しないことがあります。
実数データに対するDFTの結果の解釈
- トラブルシューティング
- 振幅スペクトルを計算する際には、複素数の絶対値 (
numpy.abs()
) を使用します。 - 位相スペクトルを計算する際には、複素数の偏角 (
numpy.angle()
) を使用します。 - 通常、実数信号の周波数解析では、正の周波数成分(配列の最初の半分程度)のみを考慮することが多いです。
numpy.fft.fftfreq()
で得られる周波数軸と対応させて解釈します。
- 振幅スペクトルを計算する際には、複素数の絶対値 (
- 原因
実数信号のDFT結果は、ナイキスト周波数を中心に共役対称性(complex conjugate symmetry)を持ちます。つまり、正の周波数成分と負の周波数成分は互いに複素共役の関係にあります。 - 問題
実数データをDFTすると、結果は複素数の配列になります。その解釈に混乱することがあります。
fftshift と ifftshift の誤用
- トラブルシューティング
fftshift
を適用したのは、表示や解析のためだけであることを意識し、逆DFTを行う前には必ずifftshift
を適用して配列の順序を元に戻します。
- 原因
fftshift
は、DFTの結果の直流成分(0 Hz)を配列の中央に移動させるための関数であり、DFT自体の結果を変更するものではありません。逆変換を行う際は、元の順序に戻す必要があります。 - 問題
周波数スペクトルの表示や解析のためにnumpy.fft.fftshift()
を使用した後、逆変換を行う前にnumpy.fft.ifftshift()
で元の順序に戻さないと、正しい時間領域の信号が得られません。
サンプリング周波数と周波数軸の不一致
- トラブルシューティング
- サンプリング周波数 (
fs
) を正確に把握し、サンプリング間隔d = 1/fs
を計算します。 numpy.fft.fftfreq(n, d)
を使用して、DFTの結果に対応する正しい周波数軸を生成します。ここでn
はDFTの点数(通常は入力データの長さ)です。
- サンプリング周波数 (
- 原因
numpy.fft.fftfreq()
関数を使用して正しい周波数軸を生成する必要があります。この関数には、DFTの点数とサンプリング間隔(またはサンプリング周波数)を正しく与える必要があります。 - 問題
DFTの結果から得られる周波数成分と、元の信号のサンプリング周波数が正しく対応付けられていないと、周波数解析の結果を誤って解釈する可能性があります。
- トラブルシューティング
- 通常、これらの誤差は実用上無視できる程度であることが多いですが、非常に高精度な計算が必要な場合は注意が必要です。
- 結果を比較する際には、厳密な等価性ではなく、ある程度の許容範囲内で比較することを検討してください (
numpy.allclose()
などを使用)。
- 原因
コンピュータでの浮動小数点数の表現には限界があるため、完全に正確な計算は難しい場合があります。 - 問題
DFTや逆DFTの計算は浮動小数点演算で行われるため、わずかな誤差が生じることがあります。
例1: 単純なサイン波のDFTと振幅スペクトルの表示
この例では、ある周波数のサイン波を生成し、そのDFTを計算して振幅スペクトルを表示します。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
fs = 100 # サンプリング周波数 (Hz)
T = 1 # サンプリング時間 (秒)
t = np.linspace(0, T, int(fs * T), endpoint=False)
f0 = 5 # 信号の周波数 (Hz)
# サイン波の生成
signal = np.sin(2 * np.pi * f0 * t)
# DFTの実行
fft_result = np.fft.fft(signal)
# 周波数軸の作成
freq = np.fft.fftfreq(signal.size, d=1/fs)
# 振幅スペクトルの計算 (複素数の絶対値)
amplitude_spectrum = np.abs(fft_result)
# プロット
plt.figure(figsize=(8, 6))
plt.plot(freq[:signal.size//2], amplitude_spectrum[:signal.size//2]) # 正の周波数成分のみ表示
plt.title(f"{f0} Hz サイン波の振幅スペクトル")
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅")
plt.grid(True)
plt.show()
解説
- パラメータ設定
サンプリング周波数 (fs
)、サンプリング時間 (T
)、信号の周波数 (f0
) を設定します。 - 時間軸の作成
np.linspace()
を使って、サンプリング時間に対応する時間軸の配列t
を作成します。endpoint=False
は、最後のサンプルを含めないようにするためです。 - サイン波の生成
指定した周波数f0
のサイン波の信号を生成します。 - DFTの実行
np.fft.fft(signal)
で信号のDFTを計算します。結果は複素数の配列fft_result
に格納されます。 - 周波数軸の作成
np.fft.fftfreq(signal.size, d=1/fs)
を使って、DFTの結果に対応する周波数軸の配列freq
を作成します。signal.size
はDFTの点数、d=1/fs
はサンプリング間隔です。 - 振幅スペクトルの計算
np.abs(fft_result)
でDFT結果の各複素数の絶対値を計算し、振幅スペクトルを得ます。 - プロット
matplotlib.pyplot
を使って、周波数軸に対して振幅スペクトルをプロットします。正の周波数成分のみを表示するために、配列の最初の半分 ([:signal.size//2]
) をスライスしています。
例2: 複数の周波数成分を含む信号のDFT
この例では、異なる2つの周波数成分を持つ信号を生成し、そのDFT結果からそれぞれの周波数を検出します。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
fs = 200 # サンプリング周波数 (Hz)
T = 2 # サンプリング時間 (秒)
t = np.linspace(0, T, int(fs * T), endpoint=False)
f1 = 10 # 1つ目の信号の周波数 (Hz)
f2 = 30 # 2つ目の信号の周波数 (Hz)
# 信号の生成
signal = np.sin(2 * np.pi * f1 * t) + 0.7 * np.sin(2 * np.pi * f2 * t)
# DFTの実行
fft_result = np.fft.fft(signal)
# 周波数軸の作成
freq = np.fft.fftfreq(signal.size, d=1/fs)
# 振幅スペクトルの計算
amplitude_spectrum = np.abs(fft_result)
# プロット
plt.figure(figsize=(10, 6))
plt.plot(freq[:signal.size//2], amplitude_spectrum[:signal.size//2])
plt.title(f"{f1} Hz と {f2} Hz の合成波の振幅スペクトル")
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅")
plt.grid(True)
plt.show()
# ピーク周波数の検出
peak_indices = np.argsort(amplitude_spectrum[:signal.size//2])[::-1][:2] # 上位2つのピークのインデックス
peak_frequencies = freq[:signal.size//2][peak_indices]
print(f"検出されたピーク周波数: {peak_frequencies} Hz")
解説
- 複数のサイン波の合成
異なる周波数 (f1
,f2
) と振幅を持つ2つのサイン波を足し合わせて合成信号を作成します。 - DFTと振幅スペクトルの計算
例1と同様にDFTを実行し、振幅スペクトルを計算します。 - プロット
振幅スペクトルをプロットすると、元の信号に含まれていた2つの周波数に対応するピークが確認できます。 - ピーク周波数の検出
np.argsort()
を使って振幅スペクトルの上位2つのピークのインデックスを取得し、それに対応する周波数をfreq
配列から抽出します。これにより、信号に含まれる主要な周波数を推定できます。
例3: 逆DFTによる信号の再構成
この例では、DFTを行った後に逆DFT (np.fft.ifft()
) を実行して、元の信号を再構成します。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定 (例1と同じ)
fs = 100
T = 1
t = np.linspace(0, T, int(fs * T), endpoint=False)
f0 = 5
# サイン波の生成
signal = np.sin(2 * np.pi * f0 * t)
# DFTの実行
fft_result = np.fft.fft(signal)
# 逆DFTの実行
reconstructed_signal = np.fft.ifft(fft_result)
# プロット
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t, signal)
plt.title("元の信号")
plt.xlabel("時間 [秒]")
plt.ylabel("振幅")
plt.subplot(2, 1, 2)
plt.plot(t, np.real(reconstructed_signal)) # 逆DFTの結果は複素数なので実部を取る
plt.title("逆DFTで再構成された信号")
plt.xlabel("時間 [秒]")
plt.ylabel("振幅")
plt.grid(True)
plt.tight_layout()
plt.show()
# 元の信号と再構成された信号の比較 (誤差)
error = np.max(np.abs(signal - np.real(reconstructed_signal)))
print(f"再構成誤差の最大値: {error}")
解説
- DFTの実行
元の信号に対してDFTを実行します。 - 逆DFTの実行
np.fft.ifft(fft_result)
を使って、DFTの結果fft_result
を時間領域の信号に戻します。逆DFTの結果は一般的に複素数として得られるため、実部 (np.real()
) を取り出して元の信号と比較します。 - プロット
元の信号と逆DFTで再構成された信号を時間領域で比較表示します。理想的には、両者はほぼ一致するはずです。 - 誤差の評価
元の信号と再構成された信号の差の絶対値の最大値を計算し、再構成の精度を評価します。浮動小数点演算の精度により、わずかな誤差が生じることがあります。
例4: fftshift
を用いた中心化された周波数スペクトルの表示
この例では、np.fft.fftshift()
を使って、DFTの結果の直流成分(0 Hz)が中央に来るように周波数スペクトルを並べ替えて表示します。
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定 (例1と同じ)
fs = 100
T = 1
t = np.linspace(0, T, int(fs * T), endpoint=False)
f0 = 5
# サイン波の生成
signal = np.sin(2 * np.pi * f0 * t)
# DFTの実行
fft_result = np.fft.fft(signal)
# 周波数軸の作成 (fftshift適用後)
freq_shifted = np.fft.fftshift(np.fft.fftfreq(signal.size, d=1/fs))
# 振幅スペクトルの計算とシフト
amplitude_spectrum = np.abs(fft_result)
amplitude_spectrum_shifted = np.fft.fftshift(amplitude_spectrum)
# プロット
plt.figure(figsize=(8, 6))
plt.plot(freq_shifted, amplitude_spectrum_shifted)
plt.title(f"{f0} Hz サイン波の中心化された振幅スペクトル")
plt.xlabel("周波数 [Hz]")
plt.ylabel("振幅")
plt.grid(True)
plt.show()
- DFTの実行と振幅スペクトルの計算
例1と同様です。 - 周波数軸のシフト
np.fft.fftshift(np.fft.fftfreq(...))
を使って、周波数軸をシフトします。これにより、負の周波数成分が左側に、正の周波数成分が右側に、そして0 Hzが中央に配置されます。 - 振幅スペクトルのシフト
同様に、np.fft.fftshift(amplitude_spectrum)
で振幅スペクトルもシフトします。 - プロット
シフトされた周波数軸に対して、シフトされた振幅スペクトルをプロットします。これにより、周波数成分の分布がより直感的に理解しやすくなります。
SciPyの scipy.fft モジュール
SciPyライブラリは、NumPyを基盤とした科学技術計算のための豊富な機能を提供しており、DFT関連の機能も scipy.fft
モジュールに用意されています。NumPyの numpy.fft
と同様の機能を提供しつつ、より高度なオプションやアルゴリズムをサポートしている場合があります。
-
例
import numpy as np from scipy.fft import fft, ifft, fftfreq import matplotlib.pyplot as plt # パラメータ設定 fs = 100 T = 1 t = np.linspace(0, T, int(fs * T), endpoint=False) f0 = 5 signal = np.sin(2 * np.pi * f0 * t) # SciPyでDFTを実行 fft_result_scipy = fft(signal) freq_scipy = fftfreq(signal.size, d=1/fs) amplitude_spectrum_scipy = np.abs(fft_result_scipy) # NumPyでDFTを実行 (比較用) fft_result_numpy = np.fft.fft(signal) freq_numpy = np.fft.fftfreq(signal.size, d=1/fs) amplitude_spectrum_numpy = np.abs(fft_result_numpy) # プロット (最初の数点のみ比較) plt.figure(figsize=(10, 6)) plt.plot(freq_numpy[:20], amplitude_spectrum_numpy[:20], label="NumPy") plt.plot(freq_scipy[:20], amplitude_spectrum_scipy[:20], linestyle='--', label="SciPy") plt.xlabel("周波数 [Hz]") plt.ylabel("振幅") plt.title("NumPy vs SciPy DFT (最初の数点)") plt.legend() plt.grid(True) plt.show()
-
利点
- 最適化された実装
SciPyは、NumPyよりもさらに最適化されたFFTアルゴリズム(例えば、FFTWライブラリの利用)を使用している場合があり、大規模なデータに対してより高速な計算が期待できることがあります。 - 高度な機能
より多くのオプション(例えば、並列処理)や、実数データに対する効率的なDFT (rfft
,irfft
) など、NumPyにはない機能を提供している場合があります。
- 最適化された実装
-
scipy.fft.fft(x, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, *, plan=None)
: 1次元DFTを行います。引数はNumPyのfft
とほぼ同様ですが、並列計算のためのworkers
オプションなどが追加されています。scipy.fft.ifft(x, n=None, axis=-1, norm=None, overwrite_x=False, workers=None, *, plan=None)
: 逆DFTを行います。scipy.fft.fftfreq(n, d=1.0)
: 周波数軸を生成します。NumPyと同様の機能です。scipy.fft.fftshift(x, axes=None)
、scipy.fft.ifftshift(x, axes=None)
: 周波数スペクトルのシフトと逆シフトを行います。
手動実装 (教育目的など)
DFTのアルゴリズムを理解するために、NumPyの機能を使わずに手動でDFTを実装することも可能です。これは通常、効率的な計算のためには推奨されませんが、DFTの原理を学ぶ上で役立ちます。
-
注意点
手動実装は計算量が O(N2) となり、データサイズが大きくなると非常に遅くなります。実用的なアプリケーションでは、NumPyやSciPyの最適化された関数を使用するべきです。 -
手動実装の例
import numpy as np import matplotlib.pyplot as plt def manual_dft(x): N = len(x) X = np.zeros(N, dtype=complex) for k in range(N): for n in range(N): X[k] += x[n] * np.exp(-2j * np.pi * k * n / N) return X # パラメータ設定 fs = 100 T = 1 t = np.linspace(0, T, int(fs * T), endpoint=False) f0 = 5 signal = np.sin(2 * np.pi * f0 * t) # 手動DFTを実行 fft_result_manual = manual_dft(signal) amplitude_spectrum_manual = np.abs(fft_result_manual) freq = np.fft.fftfreq(signal.size, d=1/fs) # NumPyのDFTを実行 (比較用) fft_result_numpy = np.fft.fft(signal) amplitude_spectrum_numpy = np.abs(fft_result_numpy) # プロット plt.figure(figsize=(10, 6)) plt.plot(freq[:signal.size//2], amplitude_spectrum_numpy[:signal.size//2], label="NumPy") plt.plot(freq[:signal.size//2], amplitude_spectrum_manual[:signal.size//2], linestyle='--', label="Manual DFT") plt.xlabel("周波数 [Hz]") plt.ylabel("振幅") plt.title("NumPy DFT vs Manual DFT") plt.legend() plt.grid(True) plt.show()
-
DFTの定義
長さ N の離散信号 x[n] (n=0,1,...,N−1) のDFT X[k] (k=0,1,...,N−1) は、以下の式で定義されます。X[k]=n=0∑N−1​x[n]e−j2πkn/N
ここで、j は虚数単位、 e−j2pikn/N は回転因子(twiddle factor)と呼ばれます。
特定の分野や用途に特化したライブラリには、DFTや関連する信号処理機能が組み込まれている場合があります。例えば、オーディオ処理のための librosa
、画像処理のための OpenCV
などです。これらのライブラリは、それぞれの分野でよく使われる処理を効率的に行うための高レベルな関数を提供しており、内部でDFTを利用していることがあります。
-
例 (librosa)
import numpy as np import librosa import librosa.display import matplotlib.pyplot as plt # サンプリングレートと時間 sr = 22050 T = 5 t = np.linspace(0, T, int(sr * T), endpoint=False) f0 = 440 # A4の周波数 # サイン波の生成 y = np.sin(2 * np.pi * f0 * t) # 短時間フーリエ変換 (STFT) を実行 (内部でDFTが使われる) stft_result = librosa.stft(y) amplitude_spectrogram = np.abs(stft_result) db_spectrogram = librosa.amplitude_to_db(amplitude_spectrogram, ref=np.max) # スペクトログラムの表示 plt.figure(figsize=(10, 4)) librosa.display.specshow(db_spectrogram, sr=sr, x_axis='time', y_axis='hz') plt.colorbar(format='%+2.0f dB') plt.title('スペクトログラム') plt.tight_layout() plt.show()
この例では、
librosa.stft()
関数が内部でDFTを利用して、時間変化する信号の周波数成分(スペクトログラム)を計算しています。