Mathematicaを使って実関数のヒルベルト変換を求める
Mathematicaを使って実関数のヒルベルト変換を求める
◼
執筆者：技術部・降旗
◼
使用したバージョン：Mathematica14.2
◼
使用した主な関数：Convolve, FourierTransform, InverseFourierTransform, Sign, Abs, Arg, ArkTan, D
はじめに
はじめに
先日ヒルベルト変換についてのお問い合わせをいただいたため、Mathematicaで実行するにはどうすれば良いかを調べてみました。
本記事ではMathematicaを使って
◼
ヒルベルト変換
◼
包絡線
◼
瞬時周波数
◼
瞬時位相
を求めてみます。
ヒルベルト変換
ヒルベルト変換
今回は実信号 をチャープ信号とし、次のように定義します。
x[t]
In[]:=
Clear["Global`*"]
In[]:=
x[t_]:=Cos+
t
4
2
t
20
In[]:=
Plot[x[t],{t,0,200}]
Out[]=
x[t]
x
x
1
πt
1
π
∞
∫
-∞
x[τ]
t-τ
によって定義されるため、 Convolve を使って畳み込みを求めてみます。
In[]:=
Convolvex[τ],τ,t,PrincipalValue->True
1
Piτ
Out[]=
ConvolveCos+,,t,PrincipalValueTrue
τ
4
2
τ
400
1
π
残念ながらこの信号では求めることができませんでした。(他の信号ではこのまま求められるものもありました。)
そこで、今度はフーリエ変換を行ってから位相をずらし、逆フーリエ変換を行う方法を試してみます。
まず、 を求めます。
そこで、今度はフーリエ変換を行ってから位相をずらし、逆フーリエ変換を行う方法を試してみます。
まず、
X[t]
In[]:=
FourierTransform[x[t],t,f]
Out[]=
(5-5)Cos-50f+100+(5+5)Cos+50f+100+(5+5)Sin-50f+100+(5-5)Sin+50f+100
25
4
2
f
25
4
2
f
25
4
2
f
25
4
2
f
続いて、 [t] を求めます。(%は一つ前に実行した結果です。計算の順序によっては同じ結果が得られないため注意が必要です。必要であれば任意の変数に代入してください。)
X
In[]:=
%*ISign[f]
Out[]=
Sign[f](5-5)Cos-50f+100+(5+5)Cos+50f+100+(5+5)Sin-50f+100+(5-5)Sin+50f+100
25
4
2
f
25
4
2
f
25
4
2
f
25
4
2
f
そして、逆フーリエ変換により [t] を求めます。(この信号では実行に少し時間がかかります。)
x
In[]:=
InverseFourierTransform[%,f,t]
Out[]=
1
2
1
400
1
400
1
200
1
20
20
2
1
200
1
20
20
2
1
200
1
20
20
2
1
200
無事に [t] が求まりました。式が少し長いため、FullSimplifyを使って簡略化を試します。
x
In[]:=
FullSimplify[%]
Out[]=
-t(100+t)Erf(50+t)+Erf(50+t)
1
2
-t(100+t)
1
400
1
200
1
20
1/4
(-1)
1
20
3/4
(-1)
In[]:=
x
1
2
-t(100+t)
1
400
1
200
1
20
1/4
(-1)
1
20
3/4
(-1)
求まった関数をプロットしてみます。
In[]:=
Plotx[t],[t],{t,0,200}
x
Out[]=
x
z[t]
In[]:=
z[t_]:=x[t]+I[t]
x
In[]:=
z[t]
Out[]=
Cos++t(100+t)Erf(50+t)+Erf(50+t)
t
4
2
t
400
1
2
-t(100+t)
1
400
1
200
1
20
1/4
(-1)
1
20
3/4
(-1)
包絡線
包絡線
包絡線は の絶対値として求められるため、Absを使用します。
z[t]
In[]:=
envelope[t_]:=Abs[z[t]]
関数をプロットしてみます。
In[]:=
Plot[{x[t],envelope[t]},{t,0,200}]
Out[]=
式の形を確認してみます。
In[]:=
envelope[t]
Out[]=
AbsCos++t(100+t)Erf(50+t)+Erf(50+t)
t
4
2
t
400
1
2
-t(100+t)
1
400
1
200
1
20
1/4
(-1)
1
20
3/4
(-1)
In[]:=
FullSimplify[%]
Out[]=
1
2
1
400
1
200
1
20
1/4
(-1)
1
20
3/4
(-1)
tが実数であると仮定すると、もう少しシンプルな形になります。
FullSimplify[%,Assumptions->t∈Reals]
瞬時位相
瞬時位相
瞬時位相 θ[t] は
の式で求められるため、ArcTanを使って求めます。
関数をプロットします。
なお、ArcTanを使わずに、Argを使ってz[t]の偏角を求める方法もあります。
ただし、こちらの方法では式にArgが残ってしまいます。
Argは微分ができないため、次の瞬時周波数を求めることができません。
Argは微分ができないため、次の瞬時周波数を求めることができません。
瞬時周波数
瞬時周波数
瞬時周波数 f[t] は
の式によって求められるため、Dを使って求めます。
求まった関数をプロットします。
おわりに
おわりに
今回はヒルベルト変換を使って実信号を複素時間信号に変換し、包絡線や瞬時周波数を求めてみました。
他のソフトでは、数値データに対してヒルベルト変換を行う例は見つかりましたが、数式に対してヒルベルト変換を行う例は見当たらなかったため、この辺りはMathematicaの強みではないかと思います。
なお、いくつかの信号を試してみましたが、信号によっては途中式がかなり長くなってしまい、計算にかかる時間も長くなるケースがありました。
そのため、場合によってはFullSimplifyを使って途中で出力される式を簡約化するなどの工夫が必要になるかもしれません。
他のソフトでは、数値データに対してヒルベルト変換を行う例は見つかりましたが、数式に対してヒルベルト変換を行う例は見当たらなかったため、この辺りはMathematicaの強みではないかと思います。
なお、いくつかの信号を試してみましたが、信号によっては途中式がかなり長くなってしまい、計算にかかる時間も長くなるケースがありました。
そのため、場合によってはFullSimplifyを使って途中で出力される式を簡約化するなどの工夫が必要になるかもしれません。