#まえがき
或る日、
教「T君、測定器から出てきたデータを送っとくから自分でグラフ作ってきて」
T 「わかりました！」
某日ミーティング・・
T 「これがグラフです」
教「ここの値なんぼだった？計算した？」
(・・・？！)
教「これとおんなじ形式のデータがまだ〇個残ってるから。計算もやっといてね。」
(（ ＾ω＾）・・・オワタ)
といった経緯があり、
データ整理を自動化、高速化するためにpythonを学習し始めました。プログラミングに手を出したのはpythonが初めてでした。少しずつ勉強し、修士論文も書き上げたところで一旦、python関連の進捗を記録していこうかと思います。記事を書くのは初で、書き終わるかも不明です(3/18)。
GitHubにもサンプルデータとnotebookがあります。こちらからどうぞ。
#データを取りました。グラフ起こし頑張ろう。
T 「今からpython勉強しながらプログラム書くのでちょっと待っててもらえますか？」
教「？？ 了！」
(後に、"ちょっと"が数カ月まで膨らむ)
さて、データを開いてみます。頂いたデータはエクセルでも開くことのできるtxt形式のデータでした。セミコロン区切りなのがちょっといやらしい。以下、データの例です。
初めて触った教科書はpython入門ノートでした。この教科書に従ってAnacondaをインストールした後、spyder上でコードを書いて勉強していました。
#ファイルの読み込み
教科書にはnumpyを使ったファイルの読み書きが書かれていました。とりあえず教科通りに読み込みました。
ファイルパスはtkinterを使って通してみました。tkinterよく分かってないですが魔法の呪文のように使っていました。
import tkinter as tk
import tkinter.filedialog as fd
import numpy as np
root=tk.Tk()
root.withdraw()
path = fd.askopenfilename(
title="file---",
filetypes=[("csv","csv"),("CSV","csv")])
if path :
fileobj=np.genfromtxt(path,delimiter=";",skip_header=3)#セミコロン区切りのデータを３行飛ばして読み込む
f=fileobj[:,0]#1列目のデータ
・・・といった具合で当時は読み込んでいました。後にpandasという便利なモジュールに出会う。
PythonユーザのためのJupyter［実践］入門 を参考にnotebook環境を立ち上げました。
jupyter notebook とpandasを引っ提げて再出発。
import tkinter
from tkinter import filedialog
import pandas as pd
root = tkinter.Tk()
root.withdraw()
path = filedialog.askopenfilename(
title="file___",
filetypes=[("txt","txt"),("csv","csv")])
if path:
df = pd.read_csv(path,engine="python",header=None,sep=';',skiprows=3,index_col=0)
#可視化
pandasで読み込むとこんな感じのテーブルデータになります。
pandasのグラフ機能でサクッとグラフ化してみます。今回の実験では1列目と2列目のデータを使います。
DataFrameではインデクサ(loc,iloc)を使ってデータを加工しますが、１行もしくは１列のデータを指定すると返ってくるオブジェクトがSeriesに変わっていることには注意が必要です。
import matplotlib.pyplot as plt
df=df.iloc[:,[0,1]]
df.set_index(0,inplace=True)#インデックスを指定してdfを上書き
df.plot()
plt.show()
pandasを使ってプロットする場合、インデックスの列を自動的に横軸に取ってくれるようですので.set_inedex(列名)で予めインデックスを指定しています。
画像のような中心あたりに鋭いピークを持つ波形が現れました。端っこの方にはノイズが乗っていました。データ点数にもよりますが、ここまではエクセルでやっていました。
#scipyを使ってフィッティング
グラフ作成までは順調に進んでいましたが、計算が面倒でした。
今回の課題はピークの鋭さを評価することでした。フィッティングツールとしてはscipyのcurve_fitがググったらすぐに出てきてたので使ってみました。
先ほどのグラフにおける縦軸は入力パワーになっておりまして、単位がデシベル(dBm)で表されていました。mWに単位を直して最大値で規格化すると下図のようになりました。
df.index = df.index*pow(10,-9)
df.index.rename('freq[GHz]',inplace=True)
df['mag'] = pow(10,(df.iloc[:]-df.iloc[:].max())/10)
df['mag'].plot()
plt.show()
フィットさせたい関数は次の通りです。ローレンツ関数にベースラインを加えたものになります。x以外は全部変数です。。
$$f_{\left(x\right)}=A \left(\frac{\gamma}{\left(x-\mu\right)^2+\gamma^2}\right)+Bx+C$$
initはパラメータの初期値です。ざっくりと予測した値で
リストを作りました。
curve_fit(関数名,x,y,パラメータの初期値)で最適なパラメータoptと共分散covが求まります。
import scipy.optimize
def lorentz(x,A,mu,gamma,B,C):#フィットする関数を定義
return A*x*gamma/((x-mu)**2+gamma**2)+B*x+C
A = -np.pi*(df.index.max()-df.index.min())/20
mu = df.index.values.mean()
gamma = (df.index.max()-df.index.min())/20
B = 10
C = 0
init = [A,mu,gamma,B,C]#フィットしたいパラメータの初期値
opt, cov = scipy.optimize.curve_fit(lorentz,df.index,df['mag'],init)#フィッティング
結果をプロットします。df [列名]=〇〇で新しく列を作れるのは大変便利。Seriesオブジェクトに追加したいときには、pd.DataFrame(Seriesオブジェクト)か.reset_index()を経てDataFrame型に変換してから。
列名をそのまま凡例にしてくれるのはうれしい点。縦軸は…。。
df['fit']=lorentz(df.index,opt[0],opt[1],opt[2],opt[3],opt[4])
df.loc[:,['mag','fit']].plot()
plt.show()
いい感じにフィットできています。
求まった$\mu$と$\gamma$を使ってピークの鋭さを評価しました。
・$\mu$ : ピークの中心
・$\gamma$ : ピークの深さが半分になるときの2点の幅 の半分(:半値半幅)
になっているのでQ値(ピークの鋭さ)
$$ Q = \frac{\mu}{2\gamma} $$
を求めて課題解決です。続く。
測定データの解析②-ヒストグラムとフィッティング,lmfitのすゝめ-
#まとめ
サンプルデータはとある共振回路の周波数特性を表していました。
独学で触り始めてやっとこさ解析プログラムが完成しました。初めて動いた時はやむごとなき感動。
卒業する頃には数千個のデータを処理していました。恐ろしい。作っておいてヨカッタ…
pandasをもっともっと勉強して自在にデータを扱えるようになりたいです。