質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

87.96%

関数を波長と強度にフーリエ変換したやつを分光器のデータ(波長800nm、強度4000くらい)に最小二乗法でfittingしたいです。

受付中

回答 0

投稿

  • 評価
  • クリップ 1
  • VIEW 226

score 1

前提・実現したいこと

関数を波長と強度にフーリエ変換したやつを分光器のデータ(波長800nm、強度4000くらい)に最小二乗法でfittingしたいです。

発生している問題・エラーメッセージ

fittingがうまくいかないです。

ガウス関数が時間tに対するyの式であり、yがフーリエ変換によってできたとしたら、
それを波長に対して並べる必要があると思いますが、
どのような操作が必要かわかりません。
(問題の原因が他のところにあるかもしれません。)

該当のソースコード

#データを読み込む(xは波長、yはスペクトル、強度)とする。
#元のデータは波長800nmで強度4000ぐらいのガウス分布

from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import pandas as pd
import csv
import math

dataset = pd.read_csv('/Users/kaiseki/taishou.csv',header=None,names=['0','1','2','3','4','5','6'])
dataset.head()


fileset=[0,50,100,150,200,250] 
x = dataset['0'] #xは波長

with open('/Users/123/0.final_report/2. bunko_fitting/bunko_kaiseki.csv','w') as f:
    writer = csv.writer(f,lineterminator='\n')



    #周波数に対して同間隔で切る方法は?
    for j in range(1,2,1):

        fig = plt.figure(figsize=(15, 10), dpi=200)
        y = dataset[str(j)]
        #yは強度
        y=np.abs(y) 


        plt.subplot(1, 3, 1)

        def func(t, *params):

            #paramsの長さでフィッティングする関数の数を判別。
            num_func = int(len(params)/4)

            #ガウス関数にそれぞれのパラメータを挿入してy_listに追加。
            y_list = []
            for i in range(num_func):
                y = np.zeros_like(t)
                param_range = list(range(4*i,4*(i+1),1))
                amp = params[int(param_range[0])] #5?
                ctr = params[int(param_range[1])]
                pulsewidth = params[int(param_range[2])] #pulseのルートが本当のパルス!
                wavelength = params[int(param_range[3])] #800?


                y1 = amp*np.exp(-(t-1/2*ctr)**2/pulsewidth**2)*np.sin((2*np.pi)*(t-1/2*ctr)/wavelength) #12/25 np.pi 追加

                y0 = y1


                spectrum = np.fft.fft(y0)

                abs_spectrum = np.abs(spectrum)  # FFTの複素数結果を絶対に変換

                #abs_spectrum=abs_spectrum[0:math.floor(step/2)]


                y = y + abs_spectrum
                y_list.append(y)

            #y_listに入っているすべてのガウス関数を重ね合わせる。
            y_sum = np.zeros_like(t)
            for i in y_list:
                y_sum = y_sum + i

            #最後にバックグラウンドを追加。
            y_sum = y_sum + params[-1]

            return y_sum






        def fit_plot(t, *params):
            num_func = int(len(params)/4)
            y_list = []
            for i in range(num_func):
                y = np.zeros_like(t)
                param_range = list(range(4*i,4*(i+1),1))
                amp = params[int(param_range[0])] #5?
                ctr = params[int(param_range[1])]
                pulsewidth = params[int(param_range[2])] #pulseのルートが本当のパルス!
                wavelength = params[int(param_range[3])] #800?

                y1 = amp*np.exp(-(t-1/2*ctr)**2/pulsewidth**2)*np.sin((2*np.pi)*(t-1/2*ctr)/wavelength)#12/25 np.pi 追加

                y0 = y1 


                spectrum = np.fft.fft(y0)
                abs_spectrum = np.abs(spectrum)  # FFTの複素数結果を絶対に変換






                y = y + abs_spectrum
                y_list.append(y)
            return y_list





        #初期値のリストを作成
        #[amp,ctr,wid]
        guess = []
        if j==1:
            guess.append([50, 1e3, 3e1, 800])
        elif j==2:
            guess.append([10,0.5e3,3e4,800])
        elif j==3:
            guess.append([100,0.34e3,6e3,800])
        elif j==4:
            guess.append([100,0.5e3,6e3,800])
        elif j==5:
            guess.append([100,0.67e3,6e3,800])
        elif j==6:
            guess.append([100,0.84e3,6e3,800])




        #バックグラウンドの初期値
        background = 5

        #初期値リストの結合
        guess_total = []
        for i in guess:
            guess_total.extend(i)
        guess_total.append(background)

        popt, pcov = curve_fit(func, x, y, p0=guess_total)

        fit = func(x, *popt)

        plt.plot(x, y, "-o")
        plt.plot(x, fit , ls='-', c='black', lw=1)

        y_list = fit_plot(x, *popt)
        baseline = np.zeros_like(x) + popt[-1]
        for n,i in enumerate(y_list):
            plt.fill_between(x, i, baseline, facecolor=cm.rainbow(n/len(y_list)), alpha=0.3)

        plt.xlim(760,840)
        plt.ylim(-100,4000)

試したこと

ここに問題に対して試したことを記載してください。

補足情報(FW/ツールのバージョンなど)

ここにより詳細な情報を記載してください。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

まだ回答がついていません

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 87.96%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る

  • トップ
  • Pythonに関する質問
  • 関数を波長と強度にフーリエ変換したやつを分光器のデータ(波長800nm、強度4000くらい)に最小二乗法でfittingしたいです。