【NumPy】高速フーリエ変換 (FFT)で振幅スペクトルを計算

この記事では、Python言語とNumPyを用いて、信号データを高速フーリエ変換(FFT)して周波数スペクトルを求める方法をソースコード付きで解説します。

高速フーリエ変換

高速フーリエ変換(Fast Fourier Transform:FFT)とは、フーリエ変換を高速化したものです。
フーリエ変換とは、デジタル信号を周波数解析するのに用いる処理です。
PythonモジュールNumpyでは「numpy.fft.fft」を用いることで高速フーリエ変換を実装できます。

フーリエ変換の原理に関してはこちら
フーリエ変換の原理・意味
逆フーリエ変換の原理・意味
離散フーリエ変換の原理・意味
高速フーリエ変換の原理・意味

メソッド

記述例 説明
F = numpy.fft.fft(f) 1次元の高速フーリエ変換により、1次元配列fのデータを時間領域から周波数領域に変換します。

ソースコード

サンプルプログラムのソースコードです。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt


def main():
    # データのパラメータ
    N = 256            # サンプル数
    dt = 0.01          # サンプリング間隔
    f1, f2 = 10, 20    # 周波数
    t = np.arange(0, N*dt, dt) # 時間軸
    freq = np.linspace(0, 1.0/dt, N) # 周波数軸

    # 信号を生成(周波数10の正弦波+周波数20の正弦波+ランダムノイズ)
    f = np.sin(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + 0.3 * np.random.randn(N)

    # 高速フーリエ変換
    F = np.fft.fft(f)

    # 振幅スペクトルを計算
    Amp = np.abs(F)

    # グラフ表示
    plt.figure()
    plt.rcParams['font.family'] = 'Times New Roman'
    plt.rcParams['font.size'] = 17
    plt.subplot(121)
    plt.plot(t, f, label='f(n)')
    plt.xlabel("Time", fontsize=20)
    plt.ylabel("Signal", fontsize=20)
    plt.grid()
    leg = plt.legend(loc=1, fontsize=25)
    leg.get_frame().set_alpha(1)
    plt.subplot(122)
    plt.plot(freq, Amp, label='|F(k)|')
    plt.xlabel('Frequency', fontsize=20)
    plt.ylabel('Amplitude', fontsize=20)
    plt.grid()
    leg = plt.legend(loc=1, fontsize=25)
    leg.get_frame().set_alpha(1)
    plt.show()


if __name__ == "__main__":
    main()

実行結果

サンプルプログラム実行結果です。
■入力信号(左)と振幅スペクトル(右)

振幅スペクトルを見ると、周波数10, 20のスペクトルが大きいことがわかります。
これは、入力信号が2つの正弦波(周波数10と20)の合成波であるためです。

補足

ナイキスト周波数f_nはサンプリング定理からf_sの半分の50となります。
(サンプリング時間(dt)が0.01なのでサンプリング周波数f_sは100)

ナイキスト周波数f_nは高速フーリエ変換の結果が有効な上限周波数となります。
先程の振幅スペクトルの場合、周波数が0~50までの部分が有効な結果となります。
(周波数が0~50以外の領域の結果に関しては無視)

以上のことから、高速フーリエ変換で周波数解析を正確に行うためには、入力信号の最大周波数の2倍以上の標本化周波数でサンプリングする必要があります。

おすすめ記事

【NumPy】高速フーリエ変換でパワースペクトル解析
Python入門 サンプル集
NumPy入門 サンプル集

シェア&フォローお願いします!