ウィグナーヴィレ分布の導出について計算メモ

#はじめに
時間周波数解析の一種である**ウィグナーヴィレ分布**のメリットデメリットやスペクトログラムとの比較については[過去の記事](
https://qiita.com/dokagui_tairan/items/1e3a1ff33c72e4bff477)で述べた．
今回はその導出方法をメモしておくことにする．

勉強のために記事を書いたものであり，正しくない説明もあるかもしれないので，コメントなどいただけると助かります．



#入力信号について
取り扱う信号は解析信号$z(t)$とする．
解析信号とは負の周波数を持たないような信号であり，`matlab`では関数`hilbert`で求められる．

[hilbert - MATLAB & Simulink - MathWorks 日本](https://jp.mathworks.com/help/signal/ref/hilbert.html)

```matlab:
Fs = 1000;
t = 0:1/Fs:1;
f = 100;
s = sin(2*pi*f*t); % 実数信号
z = hilbert(s); % 解析信号
```
`matlab`ノルマを達成したので`matlab`タグつけても許されますかね・・・？

ちなみに，$s(t)$の解析信号$z(t)$は

```math:
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$であるとすると，そのスペクトルは以下のようになる．
<img width="300" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/384b6590-6617-902f-7b01-c00b33a19904.png">


このスペクトルをデルタ関数$\delta ( \cdot )$を用いて，以下のように表す．

```math:
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$となり，ぴょこっと立つ感じである．


# 瞬時スペクトル
次に，スペクトルが時間変化する場合に拡張する．スペクトルが時々刻々と変化するので，これを瞬時スペクトルと言ったりする．

<img width="300" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/cf9d216d-17e3-fea5-f41e-66fcc3da8346.png">

瞬時周波数を$f_i(t)$として，瞬時スペクトルは，

```math:
\rho _z (t,f)=\delta ( f - f_i(t) )
```

となる．これの逆フーリエ変換を求めると，

```math:
\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)$として，それを微分したものが瞬時周波数になることは[過去の記事](https://qiita.com/dokagui_tairan/items/67728c042f3cb5ef10bb)で解説した．

```math:
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)$を用いて以下のように書ける．

```math:
\phi ' (t) = \lim_{\tau \to 0} \frac{\phi (t + \frac{\tau}{2} ) - \phi (t - \frac{\tau}{2} )} {\tau}
```


図にするとこんな感じ．
<img width="300" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/6fd72c77-a81a-a838-90e9-4eb99a2c7957.png">

さて，ここで**思い切って**以下のように近似してみる．

```math:
\phi ' (t) =  \frac{\phi (t + \frac{\tau}{2} ) - \phi (t - \frac{\tau}{2} )} {\tau}
```

この近似は$\phi'(t)$が一次関数の時，つまり，$\phi(t)$が二次関数または一次関数の場合に正しくなるようである．

これを$K_z(t,\tau)$に戻して，

```math:
\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$でフーリエ変換する．

```math:
\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)$がウィグナーヴィレ分布となる．

# まとめ
今回は時間周波数表現の答えであるところの

```math:
\rho _z (t,f)=\delta ( f - f_i(t) )
```

から逆向きにウィグナーヴィレ分布を求めてみた．


ウィグナーヴィレ分布では時間みたいな変数$\tau$が出てきたが，求め方によっては周波数みたいな変数$\nu$も出てくる．
$K_z(t,\tau)$を$t \to \nu$でフーリエ変換すると曖昧度関数$A(\nu,\tau)$が求められる．曖昧度関数，$(\nu,\tau)$平面で処理をすると何か嬉しいことがあるようだが，こちらはまだ勉強中であり，後日記事を書きたいところである．
また、`matlab`実装もゆくゆくはやっていきたいが・・・？