ロゴ

プログラミング初心者がアプリ開発を目指すブログ

Scipyで実験データのノイズを除去してみる

タイトルの通り、pythonの科学技術計算モジュールを用いて、
実験データのノイズ除去に挑戦してみました。

信号解析は全くの専門外で詳しいことはわかりませんので、やり方のみ簡単に載せます。

10Hzでサンプリングしたノイズを含んだ実験データがあります。

このデータのノイズをなるべく少なくしてみます。

まずは必要なモジュールをロードします.


>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy import fftpack

numpyのndarrayで格納した実験データ(プログラム内ではdataと表記)をフーリエ変換します。
まずはfftpack.fftfreq(data.size,d=10)でフーリエ変換の横軸を取り出します。
第1引数はデータの個数,第2引数dはサンプリング周波数(ここでは10Hz)を入れます。

fftpack.fft(data)でフーリエ変換ができます。


>>> sample_freq =fftpack.fftfreq(data.size,d=10)
>>> sig_fft=fftpack.fft(data)

sample_freqには負の値も含まれているので、
numpyのwhere関数を使って,sample_freqが正の値をとる時のインデックスを取り出します。
フーリエ変換の結果を図示してみましょう。


>>> pidxs = np.where(sample_freq > 0) #sample_freqが正の値をとる時のインデックスを取り出す.
>>> freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs] #freqsに周波数, powerに強さを入れる
>>> plt.plot(freqs, power) #図示
>>> plt.show()

今回は0.025Hz以上の部分を完全にカットしてしまいます。
フーリエ変換の0.025Hzより大きい部分の強さを0にすることで高周波のノイズを除去し、逆フーリエ変換をかけます。
逆フーリエ変換では虚数部が計算されることがあるためnumpyのreal関数をかぶせます。


>>> sig_fft[np.abs(sample_freq) > 0.025] = 0  # ノイズを消してから
>>> main_sig =np.real( fftpack.ifft(sig_fft))  # 逆フーリエ変換
>>> plt.plot(data) #元データ
[]
>>> plt.plot(main_sig)  #整形したデータ
[]
>>> plt.show()


オレンジのラインが整形したデータです。
ノイズが少なくなっています。
除去する周波数帯を変化させることでより良い結果が得られそうです。

スポンサード リンク

コメント

  • じゅん より:

    初めまして。
    >>> pidxs = np.where(sample_freq > 0) #sample_freqが正の値をとる時のインデックスを取り出す.
    >>> freqs, power = sample_freq[pidxs], np.abs(sig_fft)[pidxs] #freqsに周波数, powerに強さを入れる
    2行目のpowerに関する、sig_fftは何で定義されるのでしょうか。フーリエ変換初心者のため、基本的な質問でしたら申し訳ありません。

    • hide.me.no.car1115 より:

      じゅんさん
      コメントありがとうございます。
      powerに関する、sig_fftについてですが、私のミスで定義が抜けていました。

      プログラム5行目 sig_fft = fftpack.fft(data) とするところ
      HP内では fftpack.fft(data)となっていました.

      修正させていただきます.ご指摘ありがとうございました.

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA