Hướng dẫn voigt profile fitting python - hồ sơ voigt phù hợp với python

Tôi không thực sự chắc chắn những gì bạn đang làm hoặc cố gắng làm, nhưng đây là cách tôi sẽ làm điều đó (giả sử rằng sigma và gamma giống nhau cho tất cả các đỉnh. )

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

def cauchy(x, x0, g):
    return 1. / ( np.pi * g * ( 1 + ( ( x - x0 )/ g )**2 ) )

def gauss( x, x0, s):
    return 1./ np.sqrt(2 * np.pi * s**2 ) * np.exp( - (x-x0)**2 / ( 2 * s**2 ) )

def pseudo_voigt( x, x0, s, g, a ):
    fg = 2 * s * np.sqrt( 2 * np.log(2) )
    fl = 2 * g
    f = ( fg**5 +  2.69269 * fg**4 * fl + 2.42843 * fg**3 * fl**2 + 4.47163 * fg**2 * fl**3 + 0.07842 * fg * fl**4+ fl**5)**(1./5.)
    eta = 1.36603 * ( fl / f ) - 0.47719 * ( fl / f )**2 + 0.11116 * ( f / fl )**3
    return a * ( eta * cauchy( x, x0, f) + ( 1 - eta ) * gauss( x, x0, f ) )

def all_peaks(x, mus, amps,  s, g ):
    out = 0
    for m, a in zip( mus, amps ):
        out += pseudo_voigt( x, m, s, g, a )
    return out

def res( params, xData, yData):
    mus = params[0:5]
    amp = params[5:10]

    sig = params[-3]
    gam = params[-2]
    off = params[-1]
    yth = np.fromiter( ( abs( off ) + all_peaks( x , mus, amp, sig, gam) for x in xData ), np.float )
    diff = yth - yData
    return diff

sigma, gamma = 0.007, 0.02
offset = 0.01
muList = [ 0.5, 2.6, 4.8, 6.8,  8.9 ]
ampList = [ .135 ] * 5

data = np.loadtxt( 'calibration.txt' )
x = data[ :,0 ]
y = data[ :,1 ]

sol, err = leastsq( res, muList + ampList + [sigma , gamma, offset ], args=(x, y) )
print sol
plt.xlabel( "Voltage [V]" )
plt.ylabel( "Intensity" )

plt.plot( x,y,ls='', marker='o' )
plt.plot( x, sol[-1] + all_peaks( x, sol[0:5],sol[5:10], sol[-3], sol[-2]) )
plt.show()

mà cho

[
    4.97681822e-01 2.63788309e+00 4.74796088e+00 6.83620027e+00 8.90127524e+00 
    1.28754082e-01 1.35709531e-01 1.34679136e-01 1.35460544e-01 1.39491029e-01 
    5.61700040e-03 1.93814469e-02 9.99057213e-03
]

và biểu đồ sau

Hồ sơ dòng Voigt xảy ra trong mô hình hóa và phân tích chuyển bức xạ trong khí quyển. Đó là sự kết hợp của một hồ sơ Gaussian, $ g (x; \ sigma) $ và một hồ sơ Lorentzian, $ l (x; \ gamma) $: \ start \ int _ {-\ infty}^\ infty g (x '; \ sigma) l (x-x'; \ gamma) \, \ mathrm {d} x '\ Quad \ mathrm {WHERE} \\ g (x; \ sigma) = \ frac {1} {\ sigma \ sqrt {2 \ pi}} \ exp \ left (-\ frac {x^2} {2 \ sigma^2} \ Quad l (x; \ gamma) = \ frac {\ gamma/\ pi} {x^2 + \ gamma^2}. \ end {align*} Ở đây $ \ gamma $ là nửa chiều rộng ở nửa tối đa (HWHM) của hồ sơ Lorentzian và $ \ sigma $ là độ lệch chuẩn của hồ sơ Gaussian, liên quan đến HWHM, $ \ alpha $ $ của nó , bởi $ \ alpha = \ sigma \ sqrt {2 \ ln 2} $. Về tần số, $ \ nu $, $ x = \ nu - \ nu_0 $ trong đó $ \ nu_0 $ là trung tâm dòng.

Không có biểu mẫu đóng cho cấu hình Voigt, nhưng nó có liên quan đến phần thực của hàm faddeeva, $ w (z) $ by $$ v (x; \ sigma, \ gamma) = \ frac{\ operatorname {re}} \ re {[w (z)]}} {\ sigma \ sqrt {2 \ pi}}, \;{\ Sigma \ Sqrt {2}}.$$

Chương trình bên dưới biểu thị hồ sơ Voigt cho $ \ gamma = 0,1, \ alpha = 0,1 $ và so sánh nó với các hồ sơ Gaussian và Lorentzian tương ứng.Các phương trình trên được thực hiện trong ba hàm, G, L và V được xác định trong mã bên dưới.

import numpy as np
from scipy.special import wofz
import pylab

def G(x, alpha):
    """ Return Gaussian line shape at x with HWHM alpha """
    return np.sqrt(np.log(2) / np.pi) / alpha\
                             * np.exp(-(x / alpha)**2 * np.log(2))

def L(x, gamma):
    """ Return Lorentzian line shape at x with HWHM gamma """
    return gamma / np.pi / (x**2 + gamma**2)

def V(x, alpha, gamma):
    """
    Return the Voigt line shape at x with Lorentzian component HWHM gamma
    and Gaussian component HWHM alpha.

    """
    sigma = alpha / np.sqrt(2 * np.log(2))

    return np.real(wofz((x + 1j*gamma)/sigma/np.sqrt(2))) / sigma\
                                                           /np.sqrt(2*np.pi)

alpha, gamma = 0.1, 0.1
x = np.linspace(-0.8,0.8,1000)
pylab.plot(x, G(x, alpha), ls=':', label='Gaussian')
pylab.plot(x, L(x, gamma), ls='--', label='Lorentzian')
pylab.plot(x, V(x, alpha, gamma), label='Voigt')
pylab.xlim(-0.8,0.8)
pylab.legend()
pylab.show()