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

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

新規登録して質問してみよう
ただいま回答率
85.45%
Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

1回答

7084閲覧

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

ryokkun0331

総合スコア2

Matplotlib

MatplotlibはPythonのおよび、NumPy用のグラフ描画ライブラリです。多くの場合、IPythonと連携して使われます。

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

0クリップ

投稿2020/11/19 07:59

編集2020/11/19 08:44

前提・実現したいこと

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

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

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

イメージ説明

該当のソースコード

Python

1import numpy as np 2import math 3import matplotlib.pyplot as plt 4 5a3 = np.arange(-99.0e-12,101.0e-12,1.0e-12) #-500.0e-12~500.0e-12の時間t(要素1001個) 6A = 100 7B = 217.0e-28 8z = 30.0e+3 #伝搬距離 9T = 18016836.0e-18#パルス幅(30ps) 10dz = 10000#差分近似用Δz 11 12ft0 = A * np.exp(pow(a3[35:163]/T,2)/(-2)) 13 14plt.figure("z=30km") 15plt.title("時間領域") 16plt.xlabel("t") 17plt.ylabel("青:f(t,0)",fontname="MS Gothic") 18plt.plot(a3[35:163],ft0) 19 20plt.show() 21 22 23 24#本来のプログラム 25import numpy as np 26import math 27import matplotlib.pyplot as plt 28#範囲狭める、要素から攻める .フーリエするやつ 29#定数宣言・パラメータ 30a = np.arange(-100.0e-12,150.0e-12,1.0e-12) #500.0e-12~600.0e-12の時間t(要素1100個) 31a2 = np.arange(1,199,1) #1~1000の整数(kの一行に要素が3つある行の代入に使用) 32a3 = np.arange(-99.0e-12,101.0e-12,1.0e-12) #-500.0e-12~500.0e-12の時間t(要素1001個) 33a4 = np.arange(0,1000.0e-12,1.0e-12)#グラフ用の時間t 34k = np.eye(200) #kの土台(単位行列(1001×1001)) 35A = 100 36B = 217.0e-28 37z = 30.0e+3 #伝搬距離 38T = 18016836.0e-18#パルス幅(30ps) 39dz = 10000#差分近似用Δz 40 41#p(t,z)の作成0 42p = 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))) 43#p(t,z)の(a3サイズ)の作成 44pa3 = 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))) 45#p1a3(t+dz)差分近似用の作成 46p1a3 = 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))) 47#p2a3(t-dz)差分近似用の作成 48p2a3 = 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))) 49#dp/dzの作成 50dpdz = pow(A * T * B ,2) * z * ((2 * pow(T* a3 ,2))/(pow(T ,4)+pow(B * z ,2)) - 1) * pow(pow(T ,4) 51 + pow(B * z ,2) , -1.5) * np.exp(-1 * pow(T* a3 ,2)/(pow(T ,4) + pow(B * z ,2))) 52 53pp = (p1a3 - pa3) / dz #前進差分の作成 54pp2 = (p1a3 - p2a3) /2 /dz #中心差分の作成 55 56##f(t,0)のフーリエft0f 分散 逆フーリエft0if 57#元の式f(t,0) 58ft0 = A * np.exp(pow(a3[35:163]/T,2)/(-2)) 59ft300 = (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))) 60 61 62#フーリエ変換 63ft0f = np.fft.fft(ft0) 64 65#分散用パラメータ 66N = 128#データ数(128=2^7) 67RATE = 42.44**12 #周波数(1ビットの時間幅) 68Th = 1/RATE #ビット速度(時間幅の逆数) 69dt = Th / 32 # サンプル点間隔 70s = dt * N # 1周期時間[s] 71f = 1/s # 1周期[Hz] 72#分散式 73ome = 0 74dome = 2*np.pi*f 75ft0fb = np.zeros(N, dtype=np.complex128) #分散信号 76 77for a in range(N): 78 ft0fb[a] = ft0f[a] * np.exp(1j*B*(ome**2)*z/2) 79 ome = ome + dome 80 if a == N/2 -1: 81 ome = -ome 82 83#逆フーリエ変換 84ft0if = np.fft.ifft(ft0fb) 85ft0if = ft0if.real 86 87 88plt.figure(0) 89plt.title("f(t,0)基本",fontname="MS Gothic") 90plt.xlabel("t") 91plt.ylabel("f(t,0)",fontname="MS Gothic") 92plt.plot(a3[35:163],ft0) 93 94plt.figure(1) 95plt.title("周波数領域",fontname="MS Gothic") 96plt.xlabel("t") 97plt.ylabel("フーリエ変換後",fontname="MS Gothic") 98plt.plot(a3[35:163],ft0f) 99plt.plot(a3[35:163],ft0fb) 100 101plt.figure("z=30km") 102plt.title("時間領域") 103plt.xlabel("t") 104plt.ylabel("青:f(t,30) オレンジ:分散処理後",fontname="MS Gothic") 105#plt.plot(a3[35:163],ft0) 106plt.plot(a3[35:163],ft300) 107plt.plot(a3[35:163],ft0if) 108 109plt.show() 110 111 112#dp/dzの行列を縦一列に変形 113dpdz = dpdz[:200].reshape([200,1]) 114pp = pp[:200].reshape([200,1]) 115pp2 = pp2[:200].reshape([200,1]) 116 117 118#k(単位行列(1001×1001))の2~1000行にpを繰り返し代入 ※一行に要素が3つある行 119for i in a2: 120 k[i][i-1] = p[i+1] + p[i+2] 121 k[i][i] = -(p[i+1] + 2*p[i+2] + p[i+3]) 122 k[i][i+1] =p[i+2] + p[i+3] 123 124#kの1行目と1001行目にpを代入 ※一行に要素が2つある行 125k[199][199] = -(p[200] + 2*p[201] + p[202]) 126k[199][198] = p[200] + p[201] 127k[0][0] = -(p[1] + 2*p[2] + p[3]) 128k[0][1] = p[2] + p[3] 129 130 131kin = np.linalg.inv(k) #kの逆行列kin 132 133b = dpdz #作成とb(dpdz)の作成 134b2 = pp #b2(前進差分) 135b3 = pp2 #b3(中心差分) 136 137#各配列にどれだけ要素が入っているかの確認用(shape) 138""" 139print("a") 140print(a.shape) 141print("a2") 142print(a2.shape) 143print("a3") 144print(a3.shape) 145print(k.shape) 146print("p") 147print(p.shape) 148print("dpdz") 149print(dpdz.shape) 150print("kin") 151print(kin.shape) 152""" 153 154#各行列の確認用(shape) 155np.set_printoptions(linewidth=2000,edgeitems=7,precision=4, floatmode='maxprec') 156""" 157print("k") 158print(k) 159print("a") 160print(a) 161print("a2") 162print(a3) 163 164print("dpdz") 165print(dpdz) 166print("b") 167print(b) 168""" 169print("b") 170print(b) 171print("b2") 172print(b2) 173print("b3") 174print(b3) 175 176#3.13式よりφを計算 177iso0 = kin @ b #iso0 b=dpdz 178iso0 = iso0.reshape([200,]) 179print("iso0") 180print(iso0.shape) 181print(iso0) 182 183iso = kin @ b2#iso b=前進差分 184iso = iso.reshape([200,]) 185print("iso") 186print(iso.shape) 187print(iso) 188 189iso2 = kin @ b3#iso2 b=中心差分 190iso2 = iso2.reshape([200,]) 191print("iso2") 192print(iso2.shape) 193print(iso2) 194 195#各行列のPLOt 196""" 197plt.figure(4) 198plt.xlabel("t") 199plt.ylabel("p(t,z)") 200plt.plot(a3,p[:200]) 201 202plt.figure(5) 203plt.xlabel("t") 204plt.ylabel("bとpp") 205plt.plot(a3,dpdz) 206plt.plot(a3,pp) 207plt.plot(a3,pp2) 208 209plt.figure(6) 210plt.xlabel("t") 211plt.ylabel("φ") 212plt.plot(a3,iso0) 213plt.plot(a3,iso) 214plt.plot(a3,iso2) 215 216ig = (-B * z * pow(a3,2))/(2 * pow(T,4) + 2 * pow(B * z,2)) 217 218plt.figure(7) 219plt.xlabel("t") 220plt.ylabel("位相変化厳密",fontname="MS Gothic") 221plt.plot(a3,ig) 222 223plt.figure(8) 224plt.xlabel("t") 225plt.ylabel("位相変化厳密",fontname="MS Gothic") 226plt.plot(a3,iso0) 227plt.plot(a3,iso) 228plt.plot(a3,iso2) 229plt.plot(a3,ig) 230 231plt.show() 232""" 233 234

試したこと

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

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

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

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

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

ozwk

2020/11/19 08:08 編集

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

2020/11/19 08:08

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

2020/11/19 08:57

> 最大値がy=100であることが判明したので,y=50の時のxの値を求めたいです。 式を変形すれば正のxの値が求まります。それを2倍すれば良いかと思います。
ozwk

2020/11/19 10:01 編集

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

回答1

0

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

イメージ説明

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

物理のかぎしっぽ

投稿2020/11/21 02:04

編集2020/11/21 02:07
fkubota

総合スコア42

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

まだベストアンサーが選ばれていません

会員登録して回答してみよう

アカウントをお持ちの方は

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

ただいまの回答率
85.45%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問