【NumPy】高速フーリエ変換でパワースペクトル解析

この記事では、Python言語とNumPyを用いて、高速フーリエ変換(FFT)でパワースペクトルを計算する方法をソースコード付きで解説します。

パワースペクトル解析

パワースペクトルとは、信号の振幅と周波数の関係を示す指標です。
フーリエ変換（Fourier Transform）によりパワースペクトルを求めることができます。
今回は、PythonモジュールNumPyのnumpy.fft.fftを用いて高速フーリエ変換を行い、周波数スケールで振幅と位相をグラフ表示してみました。

書式

F = numpy.fft.fft(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)
    
    # パワースペクトルの計算（振幅スペクトルの二乗）
    Pow = Amp ** 2

    # グラフ表示
    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, Pow, 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()

実行結果

サンプルプログラム実行結果です。