ガウス型パルス(自作の波形)の半値全幅を求めたい

Question

前提・実現したいこと

ガウス型パルスの半値全幅を求めたい。
もしくは、ガウス型パルス(作成した波形)の最大値の値の半分(半値)を求め
ガウスパルス状のその半値同士の幅を求めたい

最大値がy＝100であることが判明したので,y＝50の時のxの値を求めたいです。

求めたい波形→ft0(データ数128個の配列)

該当のソースコード

import numpy as np
import math
import matplotlib.pyplot as plt

a3 = np.arange(-99.0e-12,101.0e-12,1.0e-12) #-500.0e-12～500.0e-12の時間ｔ(要素1001個)
A = 100
B = 217.0e-28
z = 30.0e+3 #伝搬距離
T = 18016836.0e-18#パルス幅（30ｐｓ）
dz = 10000#差分近似用Δｚ

ft0 = A * np.exp(pow(a3[35:163]/T,2)/(-2))

plt.figure("z=30km")
plt.title("時間領域")
plt.xlabel("t")
plt.ylabel("青:f(t,0)",fontname="MS Gothic")
plt.plot(a3[35:163],ft0)

plt.show()



#本来のプログラム
import numpy as np
import math
import matplotlib.pyplot as plt
#範囲狭める、要素から攻める .フーリエするやつ
#定数宣言・パラメータ
a = np.arange(-100.0e-12,150.0e-12,1.0e-12) #500.0e-12~600.0e-12の時間ｔ(要素1100個)
a2 = np.arange(1,199,1) #1～1000の整数(kの一行に要素が３つある行の代入に使用）
a3 = np.arange(-99.0e-12,101.0e-12,1.0e-12) #-500.0e-12～500.0e-12の時間ｔ(要素1001個)
a4 = np.arange(0,1000.0e-12,1.0e-12)#グラフ用の時間ｔ
k = np.eye(200) #kの土台(単位行列(1001×1001))
A = 100
B = 217.0e-28
z = 30.0e+3 #伝搬距離
T = 18016836.0e-18#パルス幅（30ｐｓ）
dz = 10000#差分近似用Δｚ

#p(t,z)の作成0
p = pow(A * T ,2) / np.sqrt(pow(T ,4)+pow(B * z ,2)) * np.exp((-1 * pow(T* a ,2))/(pow(T ,4)+pow(B * z ,2)))
#p(t,z)の(a3サイズ)の作成
pa3 = pow(A * T ,2) / np.sqrt(pow(T ,4)+pow(B * z ,2)) * np.exp((-1 * pow(T* a3 ,2))/(pow(T ,4)+pow(B * z ,2)))
#p1a3(t+dz)差分近似用の作成
p1a3 = pow(A * T ,2) / np.sqrt(pow(T ,4)+pow(B * (z + dz) ,2)) * np.exp((-1 * pow(T* a3 ,2))/(pow(T ,4)+pow(B * (z + dz) ,2)))
#p2a3(t-dz)差分近似用の作成
p2a3 = pow(A * T ,2) / np.sqrt(pow(T ,4)+pow(B * (z - dz) ,2)) * np.exp((-1 * pow(T* a3 ,2))/(pow(T ,4)+pow(B * (z - dz) ,2)))
#dp/dzの作成
dpdz = pow(A * T * B ,2) * z * ((2 * pow(T* a3 ,2))/(pow(T ,4)+pow(B * z ,2)) - 1) * pow(pow(T ,4)
		+ pow(B * z ,2) , -1.5)  * np.exp(-1 * pow(T* a3 ,2)/(pow(T ,4) + pow(B * z ,2)))

pp = (p1a3 - pa3) / dz #前進差分の作成
pp2 = (p1a3 - p2a3) /2 /dz #中心差分の作成

##f(t,0)のフーリエft0f　分散　逆フーリエft0if
#元の式f(t,0)
ft0 = A * np.exp(pow(a3[35:163]/T,2)/(-2))
ft300 = (A*T/np.sqrt(T**2-(1j*B*z))) * np.exp((-(T*a3[35:163])**2) / ((2 * T**4) + (2 * (B*z)**2))) * np.exp((-1j * B * z * a3[35:163]**2)/ ((2 * T**4) + (2 * (B*z)**2)))


#フーリエ変換
ft0f = np.fft.fft(ft0)

#分散用パラメータ
N = 128#データ数(128=2^7）
RATE = 42.44**12 #周波数（1ビットの時間幅）
Th = 1/RATE #ビット速度（時間幅の逆数）
dt = Th / 32 #  サンプル点間隔
s = dt * N # 1周期時間[s]
f = 1/s # 1周期[Hz]
#分散式
ome = 0
dome = 2*np.pi*f
ft0fb = np.zeros(N, dtype=np.complex128) #分散信号

for a in range(N):
	ft0fb[a] = ft0f[a] * np.exp(1j*B*(ome**2)*z/2)
	ome = ome + dome
	if a == N/2 -1:
	   ome = -ome

#逆フーリエ変換
ft0if = np.fft.ifft(ft0fb)
ft0if = ft0if.real


plt.figure(0)
plt.title("f(t,0)基本",fontname="MS Gothic")
plt.xlabel("t")
plt.ylabel("f(t,0)",fontname="MS Gothic")
plt.plot(a3[35:163],ft0)

plt.figure(1)
plt.title("周波数領域",fontname="MS Gothic")
plt.xlabel("t")
plt.ylabel("フーリエ変換後",fontname="MS Gothic")
plt.plot(a3[35:163],ft0f)
plt.plot(a3[35:163],ft0fｂ)

plt.figure("z=30km")
plt.title("時間領域")
plt.xlabel("t")
plt.ylabel("青:f(t,30) オレンジ:分散処理後",fontname="MS Gothic")
#plt.plot(a3[35:163],ft0)
plt.plot(a3[35:163],ft300)
plt.plot(a3[35:163],ft0if)

plt.show()


#dp/dzの行列を縦一列に変形
dpdz = dpdz[:200].reshape([200,1])
pp = pp[:200].reshape([200,1])
pp2 = pp2[:200].reshape([200,1])


#ｋ(単位行列(1001×1001))の2～1000行にｐを繰り返し代入　※一行に要素が３つある行
for i in a2:
	k[i][i-1] = p[i+1] + p[i+2]
	k[i][i] = -(p[i+1] + 2*p[i+2] + p[i+3])
	k[i][i+1] =p[i+2] + p[i+3]

#ｋの１行目と１００１行目にｐを代入　※一行に要素が２つある行
k[199][199] = -(p[200] + 2*p[201] + p[202])
k[199][198] = p[200] + p[201]
k[0][0] = -(p[1] + 2*p[2] + p[3])
k[0][1] = p[2] + p[3]


kin = np.linalg.inv(k) #kの逆行列kin

b = dpdz #作成とｂ(dpdz)の作成
b2 = pp #b2(前進差分)
b3 = pp2 #b3(中心差分)

#各配列にどれだけ要素が入っているかの確認用(shape)
"""
print("a")
print(a.shape)
print("a2")
print(a2.shape)
print("a3")
print(a3.shape)
print(k.shape)
print("p")
print(p.shape)
print("dpdz")
print(dpdz.shape)
print("kin")
print(kin.shape)
"""

#各行列の確認用(shape)
np.set_printoptions(linewidth=2000,edgeitems=7,precision=4, floatmode='maxprec')
"""
print("k")
print(k)
print("a")
print(a)
print("a2")
print(a3)

print("dpdz")
print(dpdz)
print("b")
print(b)
"""
print("b")
print(b)
print("b2")
print(b2)
print("b3")
print(b3)

#3.13式よりφを計算
iso0 = kin @ b #iso0 ｂ＝dpdz
iso0 = iso0.reshape([200,])
print("iso0")
print(iso0.shape)
print(iso0)

iso = kin @ b2#iso　ｂ＝前進差分
iso = iso.reshape([200,])
print("iso")
print(iso.shape)
print(iso)

iso2 = kin @ b3#iso２　b＝中心差分
iso2 = iso2.reshape([200,])
print("iso2")
print(iso2.shape)
print(iso2)

#各行列のPLOt
"""
plt.figure(4)
plt.xlabel("t")
plt.ylabel("p(t,z)")
plt.plot(a3,p[:200])

plt.figure(5)
plt.xlabel("t")
plt.ylabel("bとpp")
plt.plot(a3,dpdz)
plt.plot(a3,pp)
plt.plot(a3,pp2)

plt.figure(6)
plt.xlabel("t")
plt.ylabel("φ")
plt.plot(a3,iso0)
plt.plot(a3,iso)
plt.plot(a3,iso2)

ig = (-B * z * pow(a3,2))/(2 * pow(T,4) + 2 * pow(B * z,2))

plt.figure(7)
plt.xlabel("t")
plt.ylabel("位相変化厳密",fontname="MS Gothic")
plt.plot(a3,ig)

plt.figure(8)
plt.xlabel("t")
plt.ylabel("位相変化厳密",fontname="MS Gothic")
plt.plot(a3,iso0)
plt.plot(a3,iso)
plt.plot(a3,iso2)
plt.plot(a3,ig)

plt.show()
"""

試したこと

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

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

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

Comments

[ozwk]

やりたいことと何かしらのソースコードはありますが、
何を聞きたいのかが書いていません。
質問はなんですか?

[meg_]

最大値を求めるところから一つずつ調べてみてはどうでしょうか？

[meg_]

> 最大値がy＝100であることが判明したので,y＝50の時のxの値を求めたいです。

式を変形すれば正のxの値が求まります。それを2倍すれば良いかと思います。

[ozwk]

で、質問はなんですか？
求めたいからこう書いたらエラーが出る、なぜか？
とか
求めたいけどプログラム以前にどうしたらいいかわからないから教えてほしい
とか

Answer 1

ガウス分布であれば、半値全幅(FWHM)は標準偏差で書き表すことができます。

標準偏差をnumpyなどを使って(np.std)計算し、↑の式に代入すれば簡単に求めることができます。
この関係式の導出方法は、以下のサイトなどを参考にするといいと思います。

物理のかぎしっぽ