Matplotlib Python データ解析 モモノキ＆ナノネと学習 高速フーリエ変換FFT

Pythonで高速フーリエ変換（FFT）の練習-4 フィルタリングでノイズを除去する

Pythonで高速フリーエ変換（FFT）を行う方法をモモノキ＆ナノネと一緒に学習していきます。

モモノキ＆ナノネと一緒にPythonでFFTの使い方を覚えよう（4）

FFTとIFFTを使って信号のノイズ成分を除去してみよう

前回に続き、Pythonで高速フリー変換の練習をするよ。今回はFFTとIFFTを利用して信号のノイズ除去に挑戦しよう。

練習したFFTとIFFTを使うと信号からノイズが除去できるの？やってみる。

うん、簡単にできちゃうよ。ナノネ、まずはテスト用にノイズを含んだ信号データを作ってみて。信号は正弦波に適当な乱数を足しあわせてみて。

In [1]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

# 簡単な信号の作成(正弦波 + ノイズ)
N = 128 # サンプル数
dt = 0.01 # サンプリング周期(sec)
freq = 4 # 周波数(Hz)
amp = 1 # 振幅

t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi*freq*t) + np.random.randn(N)*0.3 # 信号

# グラフ表示
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal', fontsize=14)
plt.plot(t, f)

Out[1]:

[<matplotlib.lines.Line2D at 0xad2a2f8c>]

ノイズ入りのギザギザ信号ができた。周波数4Hzの正弦波にノイズとして乱数を少し足してみた。

ノイズ入り信号の準備はOKだね。この信号からノイズ成分を除去してみよう。最初は普通にFFTを実行してみて。

In [2]:

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

In [3]:

# FFTの複素数結果を絶対値に変換
F_abs = np.abs(F)
# 振幅をもとの信号に揃える
F_abs_amp = F_abs / N * 2 # 交流成分はデータ数で割って2倍
F_abs_amp[0] = F_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸　linspace(開始,終了,分割数)

# グラフ表示（FFT解析結果）
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F_abs_amp)

Out[3]:

[<matplotlib.lines.Line2D at 0xad15292c>]

以前の回で練習したみたいに、普通にFFTを実行したけど。ノイズ入れたから小さいピークが一杯ある。

前回の復習だけど、このまま普通にIFFT（逆変換）すると元の信号に戻るよね。

In [4]:

# そのまま普通にIFFTで逆変換した場合
F_ifft = np.fft.ifft(F) # IFFT
F_ifft_real = F_ifft.real # 実数部
plt.plot(t, F_ifft_real, c="g") # グラフ

Out[4]:

[<matplotlib.lines.Line2D at 0xad0eed8c>]

IFFTの復習でやってみたら、普通に元の信号に戻った。ノイズも復活。

ノイズ除去はFFTの結果データに少しだけ加工を加えてから、IFFTで逆変換すると実現できるよ。

ノイズを除去するデータの加工方法はいろいろあるから、その一例を試してみよう。

まず、FFT結果を別の変数にコピーしておこう。

In [5]:

F2 = np.copy(F) # FFT結果コピー

F2という変数にコピーした。

次にノイズを除去する条件を決めるよ。ノイズを取り除くにはノイズに相当するデータ値をゼロにしちゃえばOKなんだけど、ノイズとする境界値は事前に決めておく必要があるよ。

今回は周波数10を超えるデータはノイズとみなして、フィルタリング処理でデータをゼロにしちゃおう。フィルタリングに使う周波数をカットオフとも言うよ。

カットオフ（今回10Hz）を超えるデータをゼロにすればいいのね。処理はループとかで個々に判定するの？

データがNumpy配列だから、簡単にフィルターリング処理できるよ。配列の[]の中に条件を指定して、条件に合致する場合だけゼロを代入するみたいな感じだよ（Numpyの詳しい使い方はドキュメントなどで調べてね）。

今回の場合、変数fqに周波数軸の情報を代入しているから、配列[]の中の条件は(fq > カットオフ値)とすればうまくフィルタリングできるよ。

In [6]:

# 周波数でフィルタリング処理
fc = 10 # カットオフ（周波数）
F2[(fq > fc)] = 0 # カットオフを超える周波数のデータをゼロにする（ノイズ除去）

できた。カットオフ値(fc)を10、FFT結果コピー配列(F2)のデータで、条件に合致する要素をゼロにした。

カットオフを超えるデータがゼロにできたか、グラフで確認してみて。

In [14]:

# フィルタリング処理したFFT結果の確認
# FFTの複素数結果を絶対値に変換
F2_abs = np.abs(F2)
# 振幅をもとの信号に揃える
F2_abs_amp = F2_abs / N * 2 # 交流成分はデータ数で割って2倍
F2_abs_amp[0] = F2_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

# グラフ表示（FFT解析結果）
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F2_abs_amp, c='r')

Out[14]:

[<matplotlib.lines.Line2D at 0xacea598c>]

うまくいったみたい。周波数10超はゼロで真っ平らになった。

続いて、カットオフでフィルタリングしたFFT結果(F2)をIFFTで逆変換してみて。

In [ ]:

# 周波数でフィルタリング（ノイズ除去）-> IFFT
F2_ifft = np.fft.ifft(F2) # IFFT

逆変換も完了。

あとの処理は最初のときと一緒だけど、IFFTの結果から実数部を取得した値は2倍しておいてね。フィルタリングでナイキスト周波数以降も全部ゼロにしちゃったから、振幅を揃えるために必要だよ。

IFFTの結果から実数部を取得して2倍した。

In [ ]:

F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す

グラフは比較しやすいように、オリジナル信号とノイズ除去信号を重ね書きでプロットしてみて。

In [15]:

# グラフ表示：オリジナルとフィルタリング（ノイズ除去）
plt.plot(t, f, label='original')
plt.plot(t, F2_ifft_real, c="r", linewidth=4, alpha=0.7, label='filtered')
plt.legend(loc='best')
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('singnal', fontsize=14)

Out[15]:

<matplotlib.text.Text at 0xacecbeec>

できた！ノイズ除去成功かも。

うまくできたね！結果的はスーム−ジング処理に近いけど、まぁいいか。

うまくいったから、もう一つ別のフィルタリングも試してみよう。

さっきは周波数でフィルタリングしたけど、次は振幅でフィルターリングを試してみよう。

In [16]:

# グラフ再表示（FFT結果・フィルタリングなし)
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.hlines(y=[0.2], xmin=0, xmax=100, colors='r', linestyles='dashed')
plt.plot(fq, F_abs_amp)

Out[16]:

[<matplotlib.lines.Line2D at 0xace94c0c>]

上のグラフを参考に、一定振幅より小さいデータはゼロにしちゃおう。たとえば赤破線より小さい振幅は全てカットとかね。

今度は振幅の値を使ってフィルターリングすればいいのね。

今だと、グラフと同じ振幅の値はF_abs_ampに代入されているよ。あとは周波数でフィルタリングしたときと同じ要領でいけそうだね。

In [10]:

F3 = np.copy(F) # FFT結果コピー

In [11]:

# 振幅強度でフィルタリング処理
F3 = np.copy(F) # FFT結果コピー
ac = 0.2 # 振幅強度の閾値
F3[(F_abs_amp < ac)] = 0 # 振幅が閾値未満はゼロにする（ノイズ除去）

In [12]:

# 振幅でフィルタリング処理した結果の確認
# FFTの複素数結果を絶対値に変換
F3_abs = np.abs(F3)
# 振幅をもとの信号に揃える
F3_abs_amp = F3_abs / N * 2 # 交流成分はデータ数で割って2倍
F3_abs_amp[0] = F3_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要

# グラフ表示（FFT解析結果）
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F3_abs_amp, c='orange')

Out[12]:

[<matplotlib.lines.Line2D at 0xacf6abcc>]

振幅の値を使ってフィルターリングできた。FFT結果は変数名F3で、振幅0.2未満のデータはカット済み。こんどはキレイサッパリだね。

グラフ右側の鏡像ピーク（ナイキスト周波数超）もカットした方がいいの？

どっちでも大丈夫だけど、カットするのも面倒だからそのままでいいよ。その代わりさっきと違って、IFFTで逆変換したあと実数部の2倍は不要だよ。

In [13]:

# 振幅強度でフィルタリング（ノイズ除去）-> IFFT
F3_ifft = np.fft.ifft(F3) # IFFT
F3_ifft_real = F3_ifft.real # 実数部の取得
# グラフ（オリジナルとフィルタリングを比較）
plt.plot(t, f, label='original')
plt.plot(t, F3_ifft_real, c="orange", linewidth=4, alpha=0.7, label='filtered')
plt.legend(loc='best')
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('singnal', fontsize=14)

Out[13]:

<matplotlib.text.Text at 0xacf8e8ac>

うまくいったみたい。乱数がカットされて、周波数4の正弦波だけの信号になった。

OKだね、今回の練習はここまで。抜粋したまとめコードを下に載せておくね。

In [18]:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(0) # 乱数seed固定

# 簡単な信号の作成(正弦波 + ノイズ)
N = 128 # サンプル数
dt = 0.01 # サンプリング周期(sec)
freq = 4 # 周波数(Hz)
amp = 1 # 振幅

t = np.arange(0, N*dt, dt) # 時間軸
f = amp * np.sin(2*np.pi*freq*t) + np.random.randn(N)*0.3 # 信号

# 高速フーリエ変換(FFT)
F = np.fft.fft(f)
F_abs = np.abs(F) # 複素数を絶対値に変換
F_abs_amp = F_abs / N * 2 # 振幅をもとの信号に揃える(交流成分2倍)
F_abs_amp[0] = F_abs_amp[0] / 2 # 振幅をもとの信号に揃える(直流成分非2倍)

# 周波数軸のデータ作成
fq = np.linspace(0, 1.0/dt, N) # 周波数軸　linspace(開始,終了,分割数)

# フィルタリング①（周波数でカット）＊＊＊＊＊＊
F2 = np.copy(F) # FFT結果コピー
fc = 10 # カットオフ（周波数）
F2[(fq > fc)] = 0 # カットオフを超える周波数のデータをゼロにする（ノイズ除去）
F2_abs = np.abs(F2) # FFTの複素数結果を絶対値に変換
F2_abs_amp = F2_abs / N * 2 # 振幅をもとの信号に揃える(交流成分2倍)
F2_abs_amp[0] = F2_abs_amp[0] / 2 # 振幅をもとの信号に揃える(直流成分非2倍)
F2_ifft = np.fft.ifft(F2) # IFFT
F2_ifft_real = F2_ifft.real * 2 # 実数部の取得、振幅を元スケールに戻す

# フィルタリング②（振幅強度でカット）＊＊＊＊＊＊
F3 = np.copy(F) # FFT結果コピー
ac = 0.2 # 振幅強度の閾値
F3[(F_abs_amp < ac)] = 0 # 振幅が閾値未満はゼロにする（ノイズ除去）
F3_abs = np.abs(F3)# 複素数を絶対値に変換
F3_abs_amp = F3_abs / N * 2 # 交流成分はデータ数で割って2倍
F3_abs_amp[0] = F3_abs_amp[0] / 2 # 直流成分（今回は扱わないけど）は2倍不要
F3_ifft = np.fft.ifft(F3) # IFFT
F3_ifft_real = F3_ifft.real # 実数部の取得

# グラフ表示
fig = plt.figure(figsize=(12, 12))

# グラフ表示
# オリジナル信号
fig.add_subplot(321) 
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal', fontsize=14)
plt.plot(t, f)

# オリジナル信号 ->FFT
fig.add_subplot(322) 
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude', fontsize=14)
plt.plot(fq, F_abs_amp)

# オリジナル信号 ->FFT ->周波数filter ->IFFT
fig.add_subplot(323) 
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal(freq.filter)', fontsize=14)
plt.plot(t, F2_ifft_real)

# オリジナル信号 ->FFT ->周波数filter
fig.add_subplot(324) 
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude(freq.filter)', fontsize=14)
# plt.vlines(x=[10], ymin=0, ymax=1, colors='r', linestyles='dashed')
plt.fill_between([10 ,100], [0, 0], [1, 1], color='g', alpha=0.2)
plt.plot(fq, F2_abs_amp)

# オリジナル信号 ->FFT ->振幅強度filter ->IFFT
fig.add_subplot(325) 
plt.xlabel('time(sec)', fontsize=14)
plt.ylabel('signal(amp.filter)', fontsize=14)
plt.plot(t, F3_ifft_real)

# オリジナル信号 ->FFT ->振幅強度filter
fig.add_subplot(326) 
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('amplitude(amp.filter)', fontsize=14)
# plt.hlines(y=[0.2], xmin=0, xmax=100, colors='r', linestyles='dashed')
plt.fill_between([0 ,100], [0, 0], [0.2, 0.2], color='g', alpha=0.2)
plt.plot(fq, F3_abs_amp)

plt.tight_layout()

次回は周波数のピーク検出をやってみよう。

ピーク検出ね、ラジャー。

次回：Pythonで高速フーリエ変換（FFT）の練習-5 周波数ピークを自動で検出

〜　モモノキ＆ナノネと学習シリーズ（高速フーリエ変換）　〜

• Pythonで高速フーリエ変換（FFT）の練習-1 簡単な信号でFFTを体験してみよう
• Pythonで高速フーリエ変換（FFT）の練習-2 信号を時間軸と周波数軸で表現する
• Pythonで高速フーリエ変換（FFT）の練習-3 逆高速フーリエ変換（IFFT）の実践
• Pythonで高速フーリエ変換（FFT）の練習-4 フィルタリングでノイズを除去する
• Pythonで高速フーリエ変換（FFT）の練習-5 周波数ピークを自動で検出

スポンサーリンク