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

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

【Python】高速フーリエ変換

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

404 NOT FOUND | Python入門速報

メソッド

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

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

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


# データのパラメータ
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()

実行結果

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

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

サンプリング間隔が0.01なのでサンプリング周波数は100Hzとなります。
高周波側(80, 90Hz)のスペクトルは、0Hz~50Hzの正の周波数(sin)に対して、-50Hz~0Hzの負の周波数(cos)に相当する成分です。
入力信号が実数の場合、低周波側と高周波側のスペクトルは対称になりますが、複素数の場合は非対称になります。

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

【FFT応用】ローパスフィルタでノイズ除去

Python + NumPyでフーリエ変換によるローパスフィルタ処理(高周波ノイズ除去)を実装してみます。
まず、周期信号(振幅1・周波数5Hzと振幅0.2・40Hzの正弦波を重ね合わせたもの)を入力し、高周波成分(振幅0.2・40Hzの正弦波)を削ってみます。

今回実装したプログラムの処理手順は以下の通りです。

処理手順
時間信号(周波数5の正弦波 + 周波数40の正弦波)を生成します。
高速フーリエ変換により、時間信号を周波数信号に変換します。
カットオフ周波数fcを超える帯域の周波数信号の値を0にします。
高速逆フーリエ変換により、周波数信号を時間信号に変換します。
処理結果をグラフに表示します。

ソースコード

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


実行結果

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

上が「処理前の時間信号・周波数信号」、下が「処理後の時間信号・周波数信号」です。
カットオフ周波数20でローパスフィルタ処理しているので、高周波成分(周波数40の正弦波)がうまく除去できていることがわかります。

【NumPy】高速フーリエ変換とローパスフィルタでノイズ除去
Python言語とNumPyを用いて、高速逆フーリエ変換とローパス処理を行い、ノイズを取り除く方法をソースコード付きで解説します。

【FFT応用】ハイパスフィルタでノイズ除去

Python + NumPyでフーリエ変換によるハイパスフィルタ処理(高周波ノイズ抽出)を実装してみます。
今回実装したプログラムの処理手順は以下の通りです。

処理手順
時間信号(周波数5の正弦波 + 周波数40の正弦波)を生成します。
高速フーリエ変換により、時間信号を周波数信号に変換します。
カットオフ周波数fc未満の帯域と、右半分の帯域の周波数信号を0にします。
高速逆フーリエ変換により、周波数信号を時間信号に変換します。
処理結果をグラフに表示します。

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

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


def main():
    # データのパラメータ
    N = 256            # サンプル数
    dt = 0.01          # サンプリング間隔
    fq1, fq2 = 5, 40    # 周波数
    fc = 20            # カットオフ周波数
    t = np.arange(0, N*dt, dt) # 時間軸
    freq = np.linspace(0, 1.0/dt, N) # 周波数軸

    # 時間信号(周波数5の正弦波 + 周波数40の正弦波)の生成
    f = np.sin(2*np.pi*fq1*t) + 0.5 * np.sin(2*np.pi*fq2*t)

    # 高速フーリエ変換(周波数信号に変換)
    F = np.fft.fft(f)

    # 正規化 + 交流成分2倍
    F = F/(N/2)
    F[0] = F[0]/2

    # 配列Fをコピー
    F2 = F.copy()

    # ハイパス処理(カットオフ周波数未満の帯域と、右半分の帯域の周波数信号を0にする)
    F2[(freq < fc)] = 0
    F2[(freq > 1/(dt*2))] = 0

    # 高速逆フーリエ変換(時間信号に戻す)
    f2 = np.fft.ifft(F2)

    # 振幅を元のスケールに戻す
    f2 = np.real(f2*N)

    # グラフ表示
    plt.figure()
    plt.rcParams['font.family'] = 'Times New Roman'
    plt.rcParams['font.size'] = 17

    # 時間信号(元)
    plt.subplot(221)
    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(222)
    plt.plot(freq, np.abs(F), 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.subplot(223)
    plt.plot(t, f2, label='f2(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(224)
    plt.plot(freq, np.abs(F2), label='|F2(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()

実行結果

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

上が「処理前の時間信号・周波数信号」、下が「処理後の時間信号・周波数信号」です。
カットオフ周波数20でハイパスフィルタ処理しているので、高周波成分(周波数40の正弦波)がうまく抽出できていることがわかります。

【NumPy】フーリエ変換とハイパスフィルタでノイズ抽出
Python言語とNumPyを用いて、高速逆フーリエ変換とローパス処理を行い、ノイズを取り出す方法をソースコード付きで解説します。

コメント

  1. fuku より:

    質問よろしいでしょうか。
    フーリエ変換を勉強中にこのサイトに行き当たりました。
    振幅スペクトルの図についてなのですが、なぜ周波数80、90にもピークが出るのでしょうか。初歩的な質問ですみませんが、ご回答いただけますと幸いです。

    • 管理人 より:

      ※fuku様
      コメントありがとうございます。
      サンプリング定理(測定波形の最大周波数の倍の周波数以上でサンプリングしなければならない)を満たさないために、
      エイリアシング(折返し雑音)という現象が発生するためです。

      • fuku より:

        管理人様

        ありがとうございます!エイリアシングについて調べてみます。
        お礼が遅れて申し訳ありません。

        • 774 より:

          10Hzの正弦波 + 20Hzの正弦波 + 雑音なので、最大周波数は20Hzです。
          上記の例ではサンプリング間隔が0.01なのでサンプリング周波数は100Hzです。
          エイリアスは発生しません。

          高周波側に出ているのは、
          0Hz~50Hzの正の周波数(sin)に対して、-50Hz~0Hzの負の周波数(cos)に相当する成分です。
          入力が実数の場合は対称になりますが、複素数の場合は非対称になり、単に無視すればいいものではないです。

          • 管理人 より:

            774様
            丁寧にご教示いただき誠にありがとうございます。
            エイリアスに関して、完全に勘違いしておりましたので修正いたしました。

  2. yuta より:

    質問したいのですが、10hzと20hzのスペクトルで差が出るのはなぜでしょうか?