測定データの解析①-scipyフィッティングの覚書-

#まえがき
或る日、
教「T君、測定器から出てきたデータを送っとくから自分でグラフ作ってきて」
T 「わかりました！」
某日ミーティング・・
T 「これがグラフです」
教「ここの値なんぼだった？計算した？」
(・・・？！)
教「これとおんなじ形式のデータがまだ〇個残ってるから。計算もやっといてね。」
(（　＾ω＾）・・・オワタ)

といった経緯があり、
データ整理を自動化、高速化するためにpythonを学習し始めました。プログラミングに手を出したのはpythonが初めてでした。少しずつ勉強し、修士論文も書き上げたところで一旦、python関連の進捗を記録していこうかと思います。記事を書くのは初で、書き終わるかも不明です(3/18)。
GitHubにもサンプルデータとnotebookがあります。[こちらから](https://github.com/T-Ando6630/examine-fitting.git)どうぞ。

#データを取りました。グラフ起こし頑張ろう。
T 「今からpython勉強しながらプログラム書くのでちょっと待っててもらえますか？」
教「？？　了！」
(後に、"ちょっと"が数カ月まで膨らむ)

さて、データを開いてみます。頂いたデータはエクセルでも開くことのできるtxt形式のデータでした。セミコロン区切りなのがちょっといやらしい。以下、データの例です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/347339/c3da5d8b-72cc-a257-b22b-b360891717a4.png)


初めて触った教科書は[python入門ノート](https://www.amazon.co.jp/%E8%A9%B3%E7%B4%B0-Python-%E5%85%A5%E9%96%80%E3%83%8E%E3%83%BC%E3%83%88-%E5%A4%A7%E9%87%8D-%E7%BE%8E%E5%B9%B8/dp/4800711673)でした。この教科書に従ってAnacondaをインストールした後、spyder上でコードを書いて勉強していました。
#ファイルの読み込み
教科書にはnumpyを使ったファイルの読み書きが書かれていました。とりあえず教科通りに読み込みました。
ファイルパスはtkinterを使って通してみました。tkinterよく分かってないですが魔法の呪文のように使っていました。

```Python:.py
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［実践］入門](https://www.amazon.co.jp/Python%E3%83%A6%E3%83%BC%E3%82%B6%E3%81%AE%E3%81%9F%E3%82%81%E3%81%AEJupyter-%E5%AE%9F%E8%B7%B5-%E5%85%A5%E9%96%80-%E6%B1%A0%E5%86%85-%E5%AD%9D%E5%95%93/dp/4774192236) を参考にnotebook環境を立ち上げました。
jupyter notebook とpandasを引っ提げて再出発。

```Python:.ipynb
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で読み込むとこんな感じのテーブルデータになります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/347339/53f5557d-747c-7a05-9e74-90d570e49de2.png)

pandasのグラフ機能でサクッとグラフ化してみます。今回の実験では1列目と2列目のデータを使います。
DataFrameではインデクサ(loc,iloc)を使ってデータを加工しますが、１行もしくは１列のデータを指定すると返ってくるオブジェクトがSeriesに変わっていることには注意が必要です。

```Python:.ipnb
import matplotlib.pyplot as plt
df=df.iloc[:,[0,1]]
df.set_index(0,inplace=True)#インデックスを指定してdfを上書き
df.plot()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/347339/2d89aaf3-9314-5137-ef14-5ea8fb36df8b.png)

pandasを使ってプロットする場合、インデックスの列を自動的に横軸に取ってくれるようですので.set_inedex(列名)で予めインデックスを指定しています。
画像のような中心あたりに鋭いピークを持つ波形が現れました。端っこの方にはノイズが乗っていました。データ点数にもよりますが、ここまではエクセルでやっていました。
#scipyを使ってフィッティング
グラフ作成までは順調に進んでいましたが、計算が面倒でした。
今回の課題は**ピークの鋭さを評価する**ことでした。フィッティングツールとしてはscipyのcurve_fitがググったらすぐに出てきてたので使ってみました。
先ほどのグラフにおける縦軸は入力パワーになっておりまして、単位がデシベル(dBm)で表されていました。mWに単位を直して最大値で規格化すると下図のようになりました。

```Python
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()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/347339/f755454a-8e0e-6216-9d57-27da7bd8d31e.png)
フィットさせたい関数は次の通りです。ローレンツ関数にベースラインを加えたものになります。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が求まります。

```Python
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)#フィッティング
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/347339/6adb0a84-54f8-1b3f-1b40-f099f8a3fdce.png)
結果をプロットします。df [列名]=〇〇で新しく列を作れるのは大変便利。Seriesオブジェクトに追加したいときには、pd.DataFrame(Seriesオブジェクト)か.reset_index()を経てDataFrame型に変換してから。
列名をそのまま凡例にしてくれるのはうれしい点。縦軸は…。。


```Python
df['fit']=lorentz(df.index,opt[0],opt[1],opt[2],opt[3],opt[4])
df.loc[:,['mag','fit']].plot()
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/347339/56071502-94fa-cfc3-cf59-76b0f8647326.png)
いい感じにフィットできています。

求まった$\mu$と$\gamma$を使ってピークの鋭さを評価しました。
・$\mu$ : ピークの中心
・$\gamma$ : ピークの深さが半分になるときの2点の幅　の半分(:半値半幅)
になっているのでQ値(ピークの鋭さ)
$$ Q = \frac{\mu}{2\gamma} $$
を求めて課題解決です。続く。

[測定データの解析②-ヒストグラムとフィッティング,lmfitのすゝめ-](https://qiita.com/T-Ando6630/items/ba6e7510891e996dcee3)
#まとめ
サンプルデータはとある共振回路の周波数特性を表していました。
独学で触り始めてやっとこさ解析プログラムが完成しました。初めて動いた時はやむごとなき感動。
卒業する頃には数千個のデータを処理していました。恐ろしい。作っておいてヨカッタ…
pandasをもっともっと勉強して自在にデータを扱えるようになりたいです。