解析エンジニアの自動化 blog

コツコツと自動化した方法を残す blog

Python で信号のノイズ除去 〜振幅処理〜



こんにちは。
仕事の自動化にやりがいと達成感を感じるガッくんです。


この記事の目次



背景・目的


以前の記事でローパスフィルターをやりました。
Python で信号のノイズ除去 〜ローパスフィルター〜 - 解析エンジニアの自動化 blog

振幅基準でも処理できるので試してみます。

上のリンクの記事と同じ波形データ使って元の簡単な波形に近づけてみたいと思います。



動作環境


Windows 7
・winpython 64bit 3.4.4



テスト用波形

テストに使用した波形データの仕様を以下に示す。

なお、テスト用波形データはリンクの記事の波形データ作成プログラムにノイズを与えられるように改造して作成しました。

Python で複雑な波形データを作る - 解析エンジニアの自動化 blog


正弦波数:1波
  サンプリング点数:1024点
  サンプリング周期:0.05秒
  正弦波式: A × sin( 2 × π × f × t )
 
正弦波式
 テスト用波形の正弦波の式を示す。
  0.869669428796 × sin( 2 × π × 0.0675797949804 × t )

Wave for Test



プログラム

ソースコード


###############################################################################
# 波形分析プログラム - Waveform analysis program -
###############################################################################
#==============================================================================
# ファイルインポート
#==============================================================================
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
 
#==============================================================================
# 波形デジタルデータの csv 入力関数
#==============================================================================
def ReadWaveCsv():
   
    #ファイル読み込み ############################################################
    wave = np.loadtxt('C:\\Users\\UserName\\Desktop\\test-wave.csv', delimiter=',')
   
    # wave データの列スライス #####################################################
    time = wave[0:, 0] # time[sec]
    ampl = wave[0:, 1] # amplitude
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(time, ampl)
    plt.grid(True)
    plt.title('Input Wave')
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Input_Wave.jpg'
    plt.savefig(file_name)
   
    return time, ampl, file_name
 
#==============================================================================
# 波形を高速フーリエ変換する関数
#==============================================================================
def fft(time, ampl):
   
    # サンプリング周期[sec]の計算 #################################################
    sampling_cycle = time[1] - time[0]
   
    # データ数カウント ############################################################
    N = len(ampl)
   
    # 高速フーリエ変換( FFT ) ####################################################
    fft_ampl = np.fft.fft(ampl)
   
    # FFT の複素数結果を絶対に変換 ###############################################
    abs_fft_amp = np.abs(fft_ampl)
   
    # 振幅をもとの信号に揃える #####################################################
    abs_fft_amp    = abs_fft_amp    / N * 2   # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2       # 直流成分
   
    # 周波数軸のデータ作成 #######################################################
    frequency = np.linspace(0, 1.0/sampling_cycle, N) # 周波数軸 linspace(開始, 終了, 分割数)
   
    # ピーク検出 ################################################################
    maximal_idx = signal.argrelmax(abs_fft_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(frequency[:int(N/2)+1], abs_fft_amp[:int(N/2)+1])
    plt.scatter(frequency[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.grid(True)
    plt.title('Fast Fourier Transform')
    plt.xlabel('freqency[Hz]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Fast_Fourier_Transform.jpg'
    plt.savefig(file_name)
   
    return frequency, fft_ampl, file_name
 
#==============================================================================
# オペレーション・アンプリファイヤー 通称:オペアンプ
# cut_off 以下の amplitude をゼロにする
#==============================================================================
def operation_amplifier(fft_time, fft_amp, cut_off):
   
    # データ数カウント ############################################################
    N = len(fft_amp)
   
    # cut_off 処理用に実波形データに変換 ##########################################
   
    # FFT の複素数結果を絶対に変換
    abs_fft_amp = np.abs(fft_amp)
   
    # 振幅をもとの信号に揃える
    abs_fft_amp    = abs_fft_amp    / N * 2 # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2     # 直流成分
   
    # cut_off 以下の amplitude をゼロにする #######################################
    fft_amp[(abs_fft_amp <= cut_off)] = 0
   
    # グラフ用に実波形データに変換 #################################################
   
    # FFT の複素数結果を絶対に変換
    abs_fft_amp = np.abs(fft_amp)
   
    # 振幅をもとの信号に揃える
    abs_fft_amp    = abs_fft_amp    / N * 2 # 交流成分
    abs_fft_amp[0] = abs_fft_amp[0] / 2     # 直流成分
   
    # ピーク検出 ################################################################
    #maximal_idx = signal.argrelmax(abs_fft_amp[:int(N/2)+1], order=1) # 極大値インデックスの取得
    maximal_idx = signal.argrelmax(abs_fft_amp, order=1) # 極大値インデックスの取得
   
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    #plt.plot(fft_time[:int(N/2)+1], abs_fft_amp[:int(N/2)+1])
    #plt.scatter(fft_time[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.plot(fft_time, abs_fft_amp)
    plt.scatter(fft_time[maximal_idx], abs_fft_amp[maximal_idx], c='red', s=25)
    plt.grid(True)
    plt.title('Operation Amplifier')
    plt.xlabel('freqency[Hz]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Operation_Amplifier.jpg'
    plt.savefig(file_name)
   
    return fft_time, fft_amp, file_name
 
#==============================================================================
# 波形を逆高速フーリエ変換する関数
#==============================================================================
def ifft(fft_time, fft_ampl):
   
    # データ数カウント ############################################################
    N = len(fft_ampl)
   
    # 逆高速フーリエ変換
    ampl = np.fft.ifft(fft_ampl)
   
    # 実数部分のみ
    ampl_real = ampl.real
   
    # 時間軸のデータ作成 ########################################################
    frequency = fft_time[-1]
    time = np.arange(0, 1.0/frequency*N, 1.0/frequency)
    # グラフ表示 ################################################################
    plt.figure(figsize=(10, 8))
    plt.plot(time, ampl_real)
    plt.grid(True)
    plt.title('Inverse Fast Fourier Transform')
    #plt.xlim([xmin, xmax])
    #plt.ylim([-1, 2.5])
    plt.xlabel('time[sec]')
    plt.ylabel('amplitude')
   
    # グラフ出力 ################################################################
    file_name = 'Inverse_Fast_Fourier_Transform.jpg'
    plt.savefig(file_name)
   
    return time, ampl_real, file_name
 
#==============================================================================
# 実用例
#==============================================================================
# 波形の読込
time, ampl, input_wave_name = ReadWaveCsv()
 
# 高速フーリエ変換
fft_time, fft_ampl, fft_name = fft(time, ampl)
 
# オペレーション・アンプリファイヤー(cut_off 以下の amplitude をゼロにする)
fft_time, fft_ampl, opeamped_fft_name = operation_amplifier(fft_time, fft_ampl, 0.05)
 
# 逆高速フーリエ変換
ret_time, ret_ampl, ifft_name = ifft(fft_time, fft_ampl)
 
# グラフ表示
plt.figure(figsize=(10, 8))
plt.plot(time, ampl)
plt.plot(ret_time, ret_ampl)
plt.grid(True)
plt.title('Noise Cleared Wave and Original Wave')
plt.xlabel('time[sec]')
plt.ylabel('amplitude')
 
# グラフ出力
file_name = 'Noise_Cleared_Wave_and_Original_Wave.jpg'
plt.savefig(file_name)



結果

このプログラムから出力されるノイズ入り波形・極値を示した高速フーリエ変換波形・ローパスしたフーリエ変換波形・逆高速フーリエ変換波形・ノイズ除去波形とノイズ入り波形の重ね合わせをそれぞれ図1 〜 5 に示します。

図1 Input Wave


図2 Fast Fotrier Transform


図3 Operation Amplifier


図4 Inverse Fast Fotrier Transform


図5 Noise Cleared Wave and Original Wave



コメント

この方法の方が微弱(振幅の小さい)な信号に対する処理はやりやすそうです。

ピークの振幅に対して n% 以下はカットするなど使い勝手が良さそうです。

ローパスフィルターも 1 次の周期( 1 番振幅の大きな周期)に対して上下 n% の周期帯から外れた周期をカットするなど、振幅特性と周波数特性についての使い方も何となく分かりました。



以上