はじめに
時間周波数解析の一種であるウィグナーヴィレ分布のメリットデメリットやスペクトログラムとの比較については過去の記事で述べた．
今回はその導出方法をメモしておくことにする．
勉強のために記事を書いたものであり，正しくない説明もあるかもしれないので，コメントなどいただけると助かります．
入力信号について
取り扱う信号は解析信号$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実装もゆくゆくはやっていきたいが・・・？