Python FFT & IFFT
FFT & IFFT Filter プログラム
ここまで 単純な FFT & IFFT プログラムを作成してきました。
FFT を理解できたところで簡単なロ－パスフィルタ(Low-pass filter:LPF、低域通過濾波器)
のプログラムを作成します。
ロ－パスフィルタとは特定の閾値よりも高い周波数信号を減衰させて遮断し、低域周波数のみを信号として
通過させるフィルタです。
高周波域にた雑音が存在するときなどにこれを除去するのにロ－パスフィルタを用います。
ノイズを含む入力波形の作成
入力波形は基本波形とノイズ波形により作成します。
- お決まりの作法
Pyton がインスト－ルされている場所を示めします。
#!/usr/bin/python3
Pyton3 のインスト－ル場所がわからない場合は
$ which python3
で知ることができます。
ファイルで使用する文字コ－ドを設定します。
# coding: utf-8
エディタも文字コードを認識できる形式で書かれることも多いです。
# -*- coding: utf-8 -*- - ライブラリを取り込み
取り込む必要があるライブラリは２つです。
import numpy as np
import matplotlib.pyplot as plt
本ライブラリのインスト－ルがまだの方は Python/pyplot インスト－ル を見てください。
matplotlib をインスト－ルするとついでに付属で numpy も入ります。 - 入力波形の作成
サンプル数 fs = 128
振幅 amp0, amp1 = 0.8, 0.15
周波数 freq0, freq1 = 2.0, 15.0
とします。
基本波とノイズは下記のように書くことができます。
t = np.linspace(0, 2*np.pi, fs)
base = amp0 * np.sin(freq0 * t)
# 基本波
nois = amp1 * np.sin(freq1 * t)
# ノイズ
このとき t は、
スタ－トは 0
ストップは 2π = 6.28
サンプル数 fs 128
とし等間隔に t を設定します。
サンプリング周波数を f = 128Hz とすると、
サンプリングに要する時間は 1 秒必要となり、
X 軸の 0～128 の単位は 1 Hz です。
すなわち、元波形における 0～128 Hz の各周波数成分は 1 Hz ごとに分解して表示されます。
そして、0～64 Hz と 64 ～128 Hz の各周波数成分は 64 Hz を境にして 左右が対称になります。 - 入力波形の表示
plt.plot(t, base) # 基本波表示
plt.plot(t, nois) # ノイズ表示
plt.show()
とすると２つの波形をグラフで見ることができます。
また、
signal = base + nois #基本波+ノイズ
plt.plot(t, signal) # 基本波+ノイズ 表示
plt.show()
とすると入力波形をグラフで見ることができます。
フーリエ変換(FFT)～ロ－パスフィルタ処理
- フーリエ変換(FFT)
入力波形 signal をフーリエ変換します。
f = np.fft.fft(signal) # フーリエ変換(FFT)
f_abs = np.abs(f)
plt.plot(f_abs[:int(fs/2)+1])
plt.show()
FFT の結果 f は複素数です。
実数部と虚数部を別々に分けて表示しても良いのですが、今回は絶対値を取り、 かつ、左右対称なので左半分のみを表示します。 - ロ－パスフィルタ処理
フーリエ変換(FFT)結果(後述)を見ると
freq0 = 2.0, freq1 = 15.0
の部分にピ－クが存在します。
freq1 = 15.0 の部分がノイズです。
ノイズを取り除くにはノイズに相当するデ－タを０に設定するだけです。
ノイズとする境界値は、一般的にはフィルタのカットオフと呼ばれます。
ここでは、このカットオフ fc を 14 にして、
すなわち FFT 処理後の値 f の14 以上を０でクリアします。
なお、今回のケ－スではもう少し小さな値３以上を０でクリアしても結果は同じです。
fc = 14 # cut off
f[fc:]=0 # low pass filter
f_abs = np.abs(f)
plt.plot(f_abs[:int(fs/2)+1])
plt.show()
フーリエ変換(FFT)結果(後述)を見ると
freq0 = 2.0
の部分のピ－クがのみが存在します。
- ロ－パス処理後
ロ－パス処理したデ－タを IFFT 処理して時間軸の波形に変換します。
ロ－パス処理でナイキスト周波数(128/2=64) 以上も含めてすべて０に設定してしまったので、 振幅を揃えるため numpy.fft.ifft() へ与える FFT の値を２倍にします。
そして、実数部のみ表示します。
# 逆フーリエ変換(IFFT)
f_ifft = np.fft.ifft(f*2)
f_real = np.real(f_ifft)
plt.plot(t, f_real)
plt.show()
ここまでで、FFT & IFFT Filter プログラムの個々に区切った説明をしてきました。
引き続き、以上の内容をまとめたファイルを表示します。