ウィグナーヴィレ分布の導出について計算メモ
はじめに
時間周波数解析の一種であるウィグナーヴィレ分布のメリットデメリットやスペクトログラムとの比較については過去の記事で述べた.
今回はその導出方法をメモしておくことにする.
勉強のために記事を書いたものであり,正しくない説明もあるかもしれないので,コメントなどいただけると助かります.
入力信号について
取り扱う信号は解析信号$z(t)$とする.
解析信号とは負の周波数を持たないような信号であり,matlab
では関数hilbert
で求められる.
hilbert - MATLAB & Simulink - MathWorks 日本
Fs = 1000;
t = 0:1/Fs:1;
f = 100;
s = sin(2*pi*f*t); % 実数信号
z = hilbert(s); % 解析信号
matlab
ノルマを達成したのでmatlab
タグつけても許されますかね・・・?
ちなみに,$s(t)$の解析信号$z(t)$は
z(t)=\mathcal{A}(s(t))=s(t)+j \mathcal{H}(s(t))
で求められ,この$\mathcal{H}(\cdot)$がヒルベルト変換であり,matlab
では$\mathcal{A}(\cdot)$がhilbert
となっているのはちょっともにょもにょするわね.
スペクトル
$z(t)$の周波数が$f_c$であるとすると,そのスペクトルは以下のようになる.
このスペクトルをデルタ関数$\delta ( \cdot )$を用いて,以下のように表す.
Z(f) = \delta ( f - f_c )
デルタ関数$\delta(x)$は$x \not = 0$で$\delta(x)=0$,$ \int ^ { \infty } _ {- \infty } \delta (x) dx = 1$となり,
今回は$ f = f_{c} $の時にデルタ関数の中身が$0$となり,ぴょこっと立つ感じである.
瞬時スペクトル
次に,スペクトルが時間変化する場合に拡張する.スペクトルが時々刻々と変化するので,これを瞬時スペクトルと言ったりする.
瞬時周波数を$f_i(t)$として,瞬時スペクトルは,
\rho _z (t,f)=\delta ( f - f_i(t) )
となる.これの逆フーリエ変換を求めると,
\begin{align}
K_z(t,\tau)&=\mathcal{F}^{-1}_{\tau \leftarrow f} ( \rho _z (t,f) )\\
&=\mathcal{F}^{-1}_{\tau \leftarrow f} ( \delta ( f - f_i(t) ))\\
&=\int \delta ( f - f_i(t) ) \mathrm{e}^{j 2 \pi f \tau} df\\
&= \mathrm{e}^ {j 2 \pi f_i(t) \tau}\\
&= \mathrm{e}^ {j \phi ' (t) \tau}
\end{align}
時間の変数$t$はすでに使われているので,代わりに$\tau$を用いたことに注意.
また,位相$\phi(t)$として,それを微分したものが瞬時周波数になることは過去の記事で解説した.
f_i(t) = \frac{1}{2 \pi } \frac{d \phi }{dt} (t)\\
\frac{d \phi }{dt} (t) = \phi ' (t) = 2 \pi f_i (t)
さて,$\phi ' (t)$は$\phi (t)$を用いて以下のように書ける.
\phi ' (t) = \lim_{\tau \to 0} \frac{\phi (t + \frac{\tau}{2} ) - \phi (t - \frac{\tau}{2} )} {\tau}
さて,ここで思い切って以下のように近似してみる.
\phi ' (t) = \frac{\phi (t + \frac{\tau}{2} ) - \phi (t - \frac{\tau}{2} )} {\tau}
この近似は$\phi'(t)$が一次関数の時,つまり,$\phi(t)$が二次関数または一次関数の場合に正しくなるようである.
これを$K_z(t,\tau)$に戻して,
\begin{align}
K_z(t,\tau)&= \mathrm{e}^ {j \phi ' (t) \tau}\\
&=\mathrm{e}^ {j ( \phi (t + \frac{\tau}{2} ) - \phi (t - \frac{\tau}{2}) )}\\
&=\mathrm{e}^ {j \phi (t + \frac{\tau}{2} )} \mathrm{e}^ {-j \phi (t - \frac{\tau}{2} )}\\
&=z\left(t+\frac{\tau}{2}\right) z^* \left(t-\frac{\tau}{2}\right)
\end{align}
これを$\tau \to f$でフーリエ変換する.
\begin{align}
W(t,f) &= \int ^{\infty} _{- \infty} K_z(t,\tau) \mathrm{e}^{-j 2 \pi f \tau} d\tau\\
&= \int ^{\infty} _{- \infty} z\left(t+\frac{\tau}{2}\right) z^* \left(t-\frac{\tau}{2}\right) \mathrm{e}^{-j 2 \pi f \tau} d\tau\\
\end{align}
と書けば,この$W(t,f)$がウィグナーヴィレ分布となる.
まとめ
今回は時間周波数表現の答えであるところの
\rho _z (t,f)=\delta ( f - f_i(t) )
から逆向きにウィグナーヴィレ分布を求めてみた.
ウィグナーヴィレ分布では時間みたいな変数$\tau$が出てきたが,求め方によっては周波数みたいな変数$\nu$も出てくる.
$K_z(t,\tau)$を$t \to \nu$でフーリエ変換すると曖昧度関数$A(\nu,\tau)$が求められる.曖昧度関数,$(\nu,\tau)$平面で処理をすると何か嬉しいことがあるようだが,こちらはまだ勉強中であり,後日記事を書きたいところである.
また、matlab
実装もゆくゆくはやっていきたいが・・・?
Author And Source
この問題について(ウィグナーヴィレ分布の導出について計算メモ), 我々は、より多くの情報をここで見つけました https://qiita.com/dokagui_tairan/items/f9978c40c3366bd77308著者帰属:元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .