Pythonによるオシロスコープ波形データ解析の秘訣（ノイズ分析編）
2021.12.20
Pythonでノイズ分析してみた！
今回は、オシロスコープで得られた波形中の周期的なノイズ成分を Python で解析する方法を紹介します。
解析の方法としては、周期成分を抽出できるフーリエ変換を活用していき、同時にフーリエ変換の結果を上手くグラフ化する方法も説明していきます。
今回も Python で生成した波形データサンプルを使ってデータ解析をしていきます。
波形データの周期ノイズ
設計した電気回路の現物評価を行なっていく上で、直面する代表的な障害にノイズの存在があります。
電気回路のノイズには様々なタイプがありますが、最もやっかいなものは回路内の電子部品が発するタイプです。
このタイプは、電子部品自体の動作周波数に合わせる形でノイズが発生し、関係のない回路に混入することでその回路の動作を妨害することがあります。
このようなノイズの発生源を特定するためには、その特性であるノイズ周波数を特定することが重要となっています。
その周波数を特定する方法として有用な方法がフーリエ変換です。
フーリエ変換とは
フーリエ変換とは、あらゆるデータ群は周期の異なる sin/cos の三角関数の重ね合わせでできていると想定し、それぞれの周期の三角関数がデータ群にどの程度含まれているかを数値化することができる手段です。
上記はかなり曲解した表現ではありますが、ここでは他のサイトであるような学術的な省き、実用面にフォーカスした説明をしていきます。
例えば、
のような単純な正弦波であればフーリエ変換を行うと下図のように１つの鋭いピークが得られます。
このピークが立った位置の X 座標を読み取ることで、正弦波の周波数を検出することができるのです。
三角関数の含まれないデータであったとしても周期的なパターンであれば同様にピークとして検出することができます。
これは、どのような波形であっても周期がわずかにずれた三角関数の重ね合わせとして扱われるためです。
フーリエ変換のための条件
ノイズの周波数を分析したい波形データに対してフーリエ変換を行えば、波形に含まれる周期パターンの周波数を検出できることはわかりました。
では、実際にどのようなデータでも分析することができるかといいますと、それには条件があります。
まず、検出するノイズの周波数について、元の波形のサンプリング間隔が十分細かくなければ高い周波数のノイズは検出することはできません。
また、元のデータのサンプリング数、つまり周期パターンの数を十分大きくデータ内に収めなければ正確な周波数の値を検出することができません。
波形データは有限かつ離散的なものですので、解析するデータを取得する段階でどのようなサンプリング間隔でどの程度の期間のものとするかが重要となってくるわけです。
フーリエ変換後のデータ
波形データのような離散的かつ、有限なデータをフーリエ変換する場合、これは離散フーリエ変換と呼ばれ、フーリエ変換後のデータのサンプリング間隔は以下のような式で表せます。
（フーリエ変換後のサンプリング間隔[Hz]）＝１÷（データのサンプル数）÷（元データのサンプリング間隔[sec]）
前述したようにデータ数が多ければ間隔が短くなり、より正確な周波数を特定できることがわかります。
本記事で Python を使って周波数分析結果をプロットしていく際には、X 軸をこの単位の整数倍とする点に注意しておきましょう。
Python によるフーリエ変換
上で述べたフーリエ変換を Python で行うために特別なモジュールは必要ではありません。
基本モジュールの numpy に高速フーリエ変換（FFT）のメソッドがあり、それを活用する形です。
FFTは通常のフーリエ変換に対して計算を省略した方式ですが、精度面から言っても問題はありません。
使い方は分析したいデータ群を配列として引き渡すだけとなっています。
先程述べたように、変換結果のサンプル単位がデータ数に応じた周波数となることには注意しましょう。
この記事の目的
本記事では、正弦波に2つの周期的なノイズが重畳した波形を対象に、フーリエ変換を活用してノイズの周期成分を特定します。
ここでの波形は、オシロスコープで取得されたデータを想定しており、今までの記事と同様にcsv形式のデータを読み出すことを行います。
手順としては、まず Python のコード上で csv 波形データを読み出し、次に読み出したデータをフーリエ変換しつつ、フーリエ変換後の X 座標データを作成します。
グラフだけではピーク値がよく分かりませんので、ピーク検出をコード上で行い、整数倍の周波数を除いた結果が得られるようにします。
準備
今回のコード実行に必要な環境は下記になります。
インストールが必要なもの
- Python の開発環境（Spyder など）
- Pandas モジュール（波形データ生成及び csv 読み出しに使用）
- math モジュール（波形データ生成に必要）
- random モジュール（波形データ生成に必要
筆者の実行環境
- Python（3.6）
- Spyder（3.2.6）
解析する波形データ
今回、解析に使用する波形データの例は下記のコードを実行することで作成することができます。
コード本文にもあるように、データのサンプリング間隔を10usとし、1000Hz の正弦波に対して、
- ノイズ1：周期600us、周波数1666.7Hz
- ノイズ2：周期1400us、周波数714.3Hz
の２つの波形が加算された合成波形が対象となるデータです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 | import pandas as pd import random import math import matplotlib.pyplot as plt import numpy as np X_SAMPLE = 6000 NOISE_SEED_OFF = 123 WAVE_SMP_TIME = 10*10**-6 #[sec] => 10us #sin_wave SIN_INT = 100 RANDOM_NOISE_AMP=0.1 AMP= 5.0 #noise_1 NOISE_1_INT = 60 NOISE_1_INT_OFF = -20 AMP_NOISE_1 = 1.4 SIGMA_NOISE_1 = 5 noise1_int = NOISE_1_INT*WAVE_SMP_TIME noise1_freq = np.round(1/noise1_int,1) #noise_2 NOISE_2_INT = 140 NOISE_2_INT_OFF = -35 AMP_NOISE_2 = 1.0 SIGMA_NOISE_2 = 8 noise2_int = NOISE_2_INT*WAVE_SMP_TIME noise2_freq = np.round(1/noise2_int,1) xx = [ i for i in range(X_SAMPLE)] analog_wave=[] random_noise=[] for i in range(len(xx)): random.seed(NOISE_SEED_OFF+i) random_noise.append(RANDOM_NOISE_AMP*random.randint(-100,100)/100) analog_wave.append(AMP* math.sin(2/SIN_INT*xx[i]*math.pi)+random_noise[i]) noise_1=[] cnt_1=0 for i in range(len(xx)): if i % NOISE_1_INT == 0: cnt_1+=1 noise_peak_time = NOISE_1_INT* cnt_1 + NOISE_1_INT_OFF noise_1.append(AMP_NOISE_1*math.exp(-((i- noise_peak_time)**2 )/(2*SIGMA_NOISE_1**2))) noise_2=[] cnt_2=0 for i in range(len(xx)): if i % NOISE_2_INT == 0: cnt_2+=1 noise_peak_time = NOISE_2_INT* cnt_2 + NOISE_2_INT_OFF noise_2.append(AMP_NOISE_2*math.exp(-((i- noise_peak_time)**2 )/(2*SIGMA_NOISE_2**2))) #合成 merge_wave =[analog_wave[i]+noise_1[i]+noise_2[i] for i in range(X_SAMPLE)] #ノイズ1のプロット plt.figure(figsize=(10,8)) plt.plot(xx,noise_1) plt.title("Noise1,Int=" + str( np.round(noise1_int*10**(6),1))+ "[us],Freq="+str(noise1_freq)+"[Hz]") plt.savefig("./analog_noise_1.jpg") plt.show() #ノイズ2のプロット plt.figure(figsize=(10,8)) plt.plot(xx,noise_2) plt.title("Noise2,Int=" + str( np.round(noise2_int*10**(6),1))+ "[us],Freq="+str(noise2_freq)+"[Hz]") plt.savefig("./analog_noise_2.jpg") plt.show() #合成波形のプロット plt.figure(figsize=(10,8)) plt.plot(xx,merge_wave) # plt.title("Merged Wave") plt.savefig("./analog_multi_noise.jpg") plt.show() wave_data= pd.DataFrame({"multi_noise":merge_wave}) wave_data.to_csv("analog_multi_noise.csv") |
コード本文
それでは、いよいよ分析を行うコードを実行していきます。
今回のフーリエ変換では、今までの記事のようなメモリを節約できる逐次処理は使うことができません。
波形データ全体を同時に使うことでしかフーリエ変換を行うことはできないからです。
⓪モジュールのインポート
最初は、使用するモジュールのインポートからです。
普段よく使う3つをインポートすれば問題ありません。
1 2 3 | import matplotlib.pyplot as plt import numpy as np import pandas as pd |
①波形データの読み出しとフーリエ変換後のX軸単位を計算
次にcsvの読み出しとフーリエ変換後の検出可能な周波数単位を確認していきます。
1 2 3 4 5 6 | WAVE_SMP_TIME = 10*10**-6 #[sec] => 10us wave_path = "./analog_multi_noise.csv" wave_df = pd.read_csv(wave_path,index_col=0) smp_count = len(wave_df.index.values) x_unit_fft = 1/WAVE_SMP_TIME/smp_count #[Hz] print(str(x_unit_fft)+"[Hz]") |
今回は csv データの呼び出しを Pandas で行っています。
これは、今回のように全データ一括の処理を行う場合には、シンプルかつ管理がしやすいためです。
読み出したデータのサンプル数を取得し、データのサンプリング間隔を定義した後は、前述した式に基づいてフーリエ変換後の X 軸の単位（周波数分解能）を算出します。
今回の例では10us ごとのデータとしているので、16.67Hz となりました。
ここで算出した値は分解能としての意味を持つので、この周波数の整数倍でしか検出ができないということになります。
本記事ではこのまま進めますが、実際に波形測定から行う場合、なるべくたくさんのサンプルデータを取ることで測定精度上昇につながることを留意しておいてください。
データ数が多ければ処理時間が長くなり、メモリも多く必要となりますが、高周波ノイズを捉えることが目的でないならば、サンプリング間隔を広くとることでも周波数分解能を上げられることも覚えておきましょう。
②フーリエ変換の実行
そしていよいよフーリエ変換を実行します。
今回は numpy の FFT メソッドを使っていき、同時にグラフ作成用のX軸データも生成します。
変換後のデータは複素数なので絶対値も忘れずにとりましょう。
1 2 3 4 5 6 | ##フーリエ変換後のX座標単位からグラフプロット用のリストを作成する fft_x = [i*x_unit_fft for i in range(smp_count)] ##FFT処理を実施する wave_fft = np.fft.fft(wave_df.iloc[:,0].values) ##処理後の値は複素数なので絶対値をとる wave_fft_abs=abs(wave_fft) |
③ピーク位置の算出
次に、ピーク位置の検出を行います。
Python のモジュールには自動でピークを検出するメソッドは存在しないので、処理は工夫する必要があります。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | peak_lst=[] peak_lst_raw=[] itr_ct = int(len(wave_fft_abs)/2) #FFT後のX軸で中心から右半分は負の周波数なので無視する PEAK_TH=40 #データに合わせて調整する ##変換後のリストの半分を対象にピーク検出を行う for i in range(itr_ct): #ピーク検出の仕様上、最初と最後のインデックスは計算不可なのでスキップ if i!=0 and i!= itr_ct: #簡易的に対象のFFT値が前後の値より大きければピークと判定する if wave_fft_abs[i] - wave_fft_abs[i-1]>PEAK_TH \ and wave_fft_abs[i] - wave_fft_abs[i+1]>PEAK_TH: peak_lst_raw.append([i,fft_x[i],wave_fft_abs[i]]) #最初のピーク値は必ず記録する if len(peak_lst)==0: peak_lst.append([i,fft_x[i],wave_fft_abs[i]]) #検出したピーク値が既に検出したピークの整数倍なら無視する（離散化の影響を考慮し、±1のずれも同等とみなす） elif not 0 in [i % peak_lst[ii][0] for ii in range(len(peak_lst))] \ and not 0 in [(i-1) % peak_lst[ii][0] for ii in range(len(peak_lst))] \ and not 0 in [(i+1) % peak_lst[ii][0] for ii in range(len(peak_lst))]: peak_lst.append([i,fft_x[i],wave_fft_abs[i]]) |
検出方法としてはデータの最初からひとつずつ見ていき、前後の値より大きいポイントをピークとしてその X 座標と値を記録していく形です。
ここで、
1 | itr_ct = int(len(wave_fft_abs)/2) #FFT後のX軸で中心から右半分は負の周波数なので無視する |
のように検出範囲をデータの半分としていますが、これはフーリエ変換の性質上、負の周波数が後半の半分に現れるため、それを検出しないようにしているためです。
また、ポイントとして、検出したい周波数の整数倍の周波数もピークとして出てきてしまうという特徴があります。
ノイズの解析においては、検出したいのは整数倍されていない周波数なので、for 文の中の処理で、既に検出した周波数の整数倍は無視するようにしています。
分解能の関係で、必ずしも整数倍で表れませんので、プラスマイナス1ポイントのズレも整数倍とみなすコードになっています。
そして、検出した結果はリストに追加します。
④FFTによる周波数分析結果をグラフにプロット
最後に分析結果をグラフ化します。
今回のデータの分解能からすると、負の周波数を除けば最大で50MHz まで表示されてしまうので、グラフの描画範囲を絞ります。
1 2 3 4 5 | FQ_SEARCH_ST = 10 #[Hz] FQ_SEARCH_ED = 8000 #[Hz] ##グラフ描画範囲を設定する st_index = np.abs(np.asarray(fft_x) - FQ_SEARCH_ST).argmin() ed_index = np.abs(np.asarray(fft_x) - FQ_SEARCH_ED).argmin() |
そしておなじみ matplotlib でグラフ作成を行います。
1 2 3 4 5 6 7 | ##グラフのプロットとピーク検出で検出した値を表示させる plt.figure(figsize=(10,8)) plt.plot(fft_x[st_index:ed_index],wave_fft_abs[st_index:ed_index],lw=1) for i in range(len(peak_lst)): plt.annotate(str(np.round(peak_lst[i][1],1))+"Hz",xy =(peak_lst[i][1],peak_lst[i][2])) plt.savefig("./fft_result.jpg") plt.show() |
annotate を使うことでプロット上にコメントを表示することができ、先ほど検出したピークのみに周波数を表示させることができます。
作成したグラフは上図の通りで、データ作成時に1000Hz の正弦波に重畳させた1666.7Hz と714.3Hz のノイズを検出できていることがわかります。
分解能の影響でわずかに周波数値が異なる結果となっていますが、この程度の差であれば問題なく識別できるはずです。
さらにピークの大きさからノイズの大きさも相対的に判別することができるため、影響度の大きさも比較することができます。
サンプルスクリプト
最後に今まで行ったコードをコピペで実行できるように全文を下記に記載します。
基本的には、元データのサンプリング間隔と表示したい周波数範囲さえ調整すればどのようなデータにも対応できるはずです。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 | #⓪モジュールのインポート import matplotlib.pyplot as plt import numpy as np import pandas as pd #設定値関連 WAVE_SMP_TIME = 10*10**-6 #[sec] => 10us FQ_SEARCH_ST = 10 #[Hz] FQ_SEARCH_ED = 8000 #[Hz] PEAK_TH=40 #データに合わせて調整する #①波形データの読み出しとフーリエ変換後のX軸単位を計算 wave_path = "./analog_multi_noise.csv" wave_df = pd.read_csv(wave_path,index_col=0) smp_count = len(wave_df.index.values) x_unit_fft = 1/WAVE_SMP_TIME/smp_count #[Hz] print(str(x_unit_fft)+"[Hz]") #②フーリエ変換の実行 ##フーリエ変換後のX座標単位からグラフプロット用のリストを作成する fft_x = [i*x_unit_fft for i in range(smp_count)] ##FFT処理を実施する wave_fft = np.fft.fft(wave_df.iloc[:,0].values) ##処理後の値は複素数なので絶対値をとる wave_fft_abs=abs(wave_fft) #③ピーク位置の算出 peak_lst=[] peak_lst_raw=[] itr_ct = int(len(wave_fft_abs)/2) #FFT後のX軸で中心から右半分は負の周波数なので無視する ##変換後のリストの半分を対象にピーク検出を行う for i in range(itr_ct): #ピーク検出の仕様上、最初と最後のインデックスは計算不可なのでスキップ if i!=0 and i!= itr_ct: #簡易的に対象のFFT値が前後の値より大きければピークと判定する if wave_fft_abs[i] - wave_fft_abs[i-1]>PEAK_TH \ and wave_fft_abs[i] - wave_fft_abs[i+1]>PEAK_TH: peak_lst_raw.append([i,fft_x[i],wave_fft_abs[i]]) #最初のピーク値は必ず記録する if len(peak_lst)==0: peak_lst.append([i,fft_x[i],wave_fft_abs[i]]) #検出したピーク値が既に検出したピークの整数倍なら無視する（離散化の影響を考慮し、±1のずれも同等とみなす） elif not 0 in [i % peak_lst[ii][0] for ii in range(len(peak_lst))] \ and not 0 in [(i-1) % peak_lst[ii][0] for ii in range(len(peak_lst))] \ and not 0 in [(i+1) % peak_lst[ii][0] for ii in range(len(peak_lst))]: peak_lst.append([i,fft_x[i],wave_fft_abs[i]]) #④FFTによる周波数分析結果をグラフにプロット ##グラフ描画範囲を設定する st_index = np.abs(np.asarray(fft_x) - FQ_SEARCH_ST).argmin() ed_index = np.abs(np.asarray(fft_x) - FQ_SEARCH_ED).argmin() ##グラフのプロットとピーク検出で検出した値を表示させる plt.figure(figsize=(10,8)) plt.plot(fft_x[st_index:ed_index],wave_fft_abs[st_index:ed_index],lw=1) for i in range(len(peak_lst)): plt.annotate(str(np.round(peak_lst[i][1],1))+"Hz",xy =(peak_lst[i][1],peak_lst[i][2])) plt.savefig("./fft_result.jpg") plt.show() |
こちらの記事もオススメ！
書いた人はこんな人
- 「好きを仕事にするエンジニア集団」の（株）ライトコードです！
ライトコードは、福岡、東京、大阪の3拠点で事業展開するIT企業です。
現在は、国内を代表する大手IT企業を取引先にもち、ITシステムの受託事業が中心。
いずれも直取引で、月間PV数1億を超えるWebサービスのシステム開発・運営、インフラの構築・運用に携わっています。
システム開発依頼・お見積もり大歓迎！
また、現在「WEBエンジニア」「モバイルエンジニア」「営業」「WEBデザイナー」「WEBディレクター」を積極採用中です！
インターンや新卒採用も行っております。
以下よりご応募をお待ちしております！
https://rightcode.co.jp/recruit
ライトコードの日常12月 1, 2023ライトコードクエスト〜東京オフィス歴史編〜- ITエンタメ10月 13, 2023Netflixの成功はレコメンドエンジン？
- ライトコードの日常8月 30, 2023退職者の最終出社日に密着してみた！
- ITエンタメ8月 3, 2023世界初の量産型ポータブルコンピュータを開発したのに倒産!?アダム・オズボーン
関連記事
採用情報
\ あの有名サービスに参画!? /
バックエンドエンジニア
\ クリエイティブの最前線 /
フロントエンドエンジニア
\ 世界を変える…！ /
Androidエンジニア
\ みんなが使うアプリを創る /
iOSエンジニア
