【MATLAB】ウィグナーヴィレ分布とスペクトログラムを比較する【時間周波数解析】

#はじめに
今回は時間周波数解析の一種である**ウィグナーヴィレ分布**について紹介し，スペクトログラムと比較しながらその利点と欠点について解説する．
詳しい計算を解説するとかなりのボリュームになるので，その辺は後日記事を書こうと思う．
今回は計算無しで行く．

#ウィグナーヴィレ分布とは
計算無しで行くと言ったがとりあえず定義だけ書く．ウィグナーヴィレ分布$W(t,f)$は，入力信号$x(t)$の解析信号$z(t)=\mathcal{A}(x(t))$とすると以下のように定義される．

```math:
W(t,f)=\int ^ \infty _{-\infty} z\left(t+\frac{\tau}{2} \right)z^*\left(t-\frac{\tau}{2} \right) \mathrm{e}^{-j2 \pi f \tau} d\tau
```

解析信号とは入力信号を負の周波数を持たないように変換をした複素信号である．
なお，実数信号をそのまま用いるとウィグナー分布$\mathcal{W}(t,f)$になる．

```math:
\mathcal{W}(t,f)=\int ^ \infty _{-\infty} s\left(t+\frac{\tau}{2} \right)s\left(t-\frac{\tau}{2} \right) \mathrm{e}^{-j2 \pi f \tau} d\tau
```

ウィグナーヴィレ分布の特徴として**かなりシャープな**結果を出してくれる．以下に例を挙げていこう．

# 例 1つのチャープ信号を分析してみる
以下の瞬時周波数を持つチャープ信号を分析してみる．

```math:
f_i(t)=150t+100
```

<img width="400" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/1de73fba-6afa-c226-88bb-3feea335785f.png">

<details><summary>チャープ信号の解説は↓</summary><div>

$t=0\rm{[s]}$で$f=100\rm{[Hz]}$であり，$t=1\rm{[s]}$ごとに$f=150\rm{[Hz]}$だけ周波数が高くなっていく信号である．


チャープ信号は，瞬時位相$\phi_i(t)=2 \pi f_i$を積分して$\sin$か$\cos$に突っ込めばよいので，

```math:
y(t)=\cos \left( 2 \pi \left( 150\frac{1}{2}t^2+100t \right) + \phi_0 \right) 
```

時間波形は，

```matlab:
Fs = 1000;
T = 2;
N = Fs*T;
t = ((1:N)-1)/Fs;
fi = 150*t + 100;
figure('Position' , [100 100 1000 600])
plot(t,fi);
ylim([0,Fs/2]);
xlabel('時間[s]');
ylabel('周波数[Hz]');
phi = 150/2*t.^2 + 100*t;
y = cos(2*pi*(phi));
figure('Position' , [100 100 1000 600])
plot(t,y);
ylim([-1 1]);
xlabel('時間[s]');
ylabel('振幅');
```
チャープ信号の時間波形は以下の通り．
<img width="400" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/9e36ffc1-45c9-b650-5fd0-3d3e92ca6cc0.png">


</div></details>

##スペクトログラムの結果
比較のため，一般的な時間周波数解析手法であるスペクトログラム（短時間フーリエ変換を絶対値の$2$乗したもの）で分析をする．

```matlab:
figure('Position', [100 100 1000 600])
Nw = 128;
Nol = round(Nw*0.95);
NFFT = Nw;
spectrogram(y, Nw, Nol, NFFT, Fs, 'yaxis');
```

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/9ea4a1d4-f44b-125b-5945-6b0a0ad66136.png">


##ウィグナーヴィレ分布の結果
今回は`Signal Processing Toolbox`の`wvd`を使用する．

[Wigner-Ville 分布と平滑化疑似 Wigner-Ville 分布 - MathWorks](https://jp.mathworks.com/help/signal/ref/wvd.html)

```matlab:
figure('Position', [100 100 1000 600])
wvd(y, Fs);
caxis([0 1500])
```

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/8902bc8a-0aff-148e-155e-ace381fdee8e.png">


スペクトログラムと比較して，**めちゃくちゃシャープに**結果が表示されていることが分かると思う．
**スゴイ！**


#例　2つのチャープ信号
次は以下の瞬時周波数を持つ2つのチャープ信号を分析する．

```math:
f_1(t)=100t+200\\
f_2(t)=100t+50
```

<img width="400" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/079c3eab-1777-1c02-77f8-08925594297a.png">

<details><summary>解説及び`matlab`実装は↓</summary><div>

```math:
y(t) = \cos \left( 2 \pi \left( 100\frac{1}{2}t^2+100t \right) + \phi_0 \right) 
 + \cos \left( 2 \pi \left( 100\frac{1}{2}t^2+50t \right) + \psi_0 \right) 
```


```matlab:
Fs = 1000;
T = 2;
N = Fs*T;
t = ((1:N)-1)/Fs;
%fi = 150*t + 100;
fi1 = 100*t + 200;
fi2 = 100*t + 50;
figure('Position' , [100 100 1000 600])
plot(t,fi1);
hold on
plot(t,fi2);
ylim([0,Fs/2]);
xlabel('時間[s]');
ylabel('周波数[Hz]');

%phi = 150/2*t.^2 + 100*t;
phi = 100/2*t.^2 + 200*t;
phi2 =  100/2*t.^2 + 50*t;
y = cos(2*pi*(phi)) + cos(2*pi*(phi2));
figure('Position' , [100 100 1000 600])
plot(t,y);
ylim([-1 1]);
xlabel('時間[s]');
ylabel('振幅');
% スペクトログラム
figure('Position' , [100 100 1000 600])
Nw = 128;
Nol = round(Nw*0.95);
NFFT = Nw;
spectrogram(y, Nw, Nol, NFFT, Fs, 'yaxis');
% ウィグナーヴィレ分布
figure('Position' , [100 100 1000 600])
wvd(y, Fs);
caxis([0 1500]);% そのままだと線が見えづらかったので表示レベル範囲をちょと修正
```

</div></details>

##スペクトログラムの結果
スペクトログラム
<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/2b2152fe-9c3b-6fa5-8a8e-81a38dcdff68.png">


##ウィグナーヴィレ分布の結果
<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/1b2a2ada-9742-97f6-601a-a73865b46563.png">

ん？

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/42b7bcee-1b4c-11b9-07e8-fbc383fb37b2.png">

ウィグナーヴィレ分布では**謎の線**が出てしまっている．


#例　2つのチャープ信号その2
気を取り直して，次はこちらの瞬時周波数を持つチャープ信号を分析しよう．

```math:
f_1(t)=100t+200\\
f_2(t)=-150t+400
```

<img width="400" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/db3cdad9-5fe5-62f4-c11e-e3d39290062a.png">



<details><summary>解説及び`matlab`実装は↓</summary><div>

```math:
y(t) = \cos \left( 2 \pi \left( 100\frac{1}{2}t^2+100t \right) + \phi_0 \right) 
 + \cos \left( 2 \pi \left( -150\frac{1}{2}t^2+400t \right) + \psi_0 \right) 
```


```matlab:
phi = 100/2*t.^2 + 200*t;
phi2 =  100/2*t.^2 + 50*t;
y = cos(2*pi*(phi)) + cos(2*pi*(phi2));
% とすれば，あとは上と同じ．
%wvdの表示範囲は↓とした．
%caxis([0 1000])
```

</div></details>

##スペクトログラムの結果

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/6739d73b-3508-6e5e-78b1-463997a99e0e.png">


##ウィグナーヴィレ分布の結果

<img width="600" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/19188478-89a7-d2b3-fb8a-c72a8a27dca4.png">


ん～～～？？？
![wigner3_2._uppng.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/0782545c-9d4d-b5c7-af60-4b908077a880.png)
$2$線が交わる部分を拡大してみた．
**縞模様**のようなものが見えることが分かる．




#ウィグナーヴィレ分布で発生する謎成分は何？
ウィグナーヴィレ分布で発生する謎のスペクトル成分（アーティファクト，偽像）．
これを**クロス項**という．

計算は後日記事を書くことにするが，$z(t)=z_1(t)+z_2(t)$のウィグナーヴィレ分布は以下のようになる．

```math:
W_z=W_{z_1}+W_{z_2}+2\rm{Re}(W_{z_1 z_2})
```

ここで

```math:
W_{z_1 z_2}=\int ^ \infty _{-\infty} z_1\left(t+\frac{\tau}{2} \right)z_2^*\left(t-\frac{\tau}{2} \right) \mathrm{e}^{-j2 \pi f \tau} d\tau
```

この$W_{z_1 z_2}$が**クロス項**である．


<img width="300" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/1b2a2ada-9742-97f6-601a-a73865b46563.png">
再掲．この図を見ると分かりやすいが，クロス項は$2$信号の**中間の位置**に発生する．



ウィグナーヴィレ分布はこのように，複数の信号があった場合に，クロス項が発生することが厄介な点である．
今回の例は単純なので分かりやすいが，実際の信号だとどれが本物の信号でどれがクロス項か判断がつかないことがある．
**界隈**ではこれを**結果の解釈に専門性が問われる**と言ったりする．


#ウィグナーヴィレ分布のその他の問題点

##計算が遅い
メモリがかなり多めに必要だし，処理時間がかなりかかるようである．

##位相の情報が失われる
逆変換が難しい．ウィグナーヴィレ分布を求めた後で時間周波数平面で雑音除去などの処理をして，それを再び音声データに戻すことが難しい．


#クロス項の問題を解決する方法はあるの？
クロス項を抑制する方法は色々と研究されている．
<img width="300" alt="" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373168/0782545c-9d4d-b5c7-af60-4b908077a880.png">
再掲．これを見ると縞模様が出来ており，言い方を変えると**クロス項は時間方向や周波数方向に振動している**わけである．
よって，これに**ローパスフィルタをかける**とクロス項を抑制することが期待できる．
$(t-f)$平面ではなく，**曖昧度関数**$(\tau-\nu)$平面で処理をすることもある．

この辺を一般化した話が**コーエンクラス**である．

結果はボケボケだけど安定しているスペクトログラムと，結果はキレキレだけどクロス項が発生しちゃうウィグナーヴィレ分布．
キレキレの結果に平滑化フィルタをかけて結果をボカすのは**元の木阿弥**なのではと思うし実際平滑化フィルタとして**ガウシアン**をかけると結果が**スペクトログラムになってしまう**らしいのだが，このあたりは勉強中なので後日記事を書きたい。