【NumPy】高速フーリエ変換とローパスフィルタでノイズ除去

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

フーリエ変換によるノイズ除去(ローパスフィルタ処理)

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.2 * 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 = 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の正弦波)がうまく除去できていることがわかります。

関連記事
1 Python入門 基本文法
2 NumPy入門 サンプル集
参考文献 Python NumPy SciPy サンプルコード: フーリエ変換処理 その 3
関連記事