MATLAB実装から逃げるな
時間周波数解析手法のウィグナーヴィレ分布については以下の記事を書いてきた．
【MATLAB】ウィグナーヴィレ分布とスペクトログラムを比較する【時間周波数解析】
ウィグナーヴィレ分布の導出について計算メモ
今回はついにmatlabと和解することにし，ウィグナーヴィレ分布をmatlabで実装してみることにする．
matlabなんも分からんし，信号処理も雰囲気でやっているため，実装についてはあまり参考にしてはいけない点に注意．
ウィグナーヴィレ分布おさらい
ウィグナーヴィレ分布は以下のように書ける．
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
Signal Kernel
ここでsignal kernel $K_z$を以下のように定義する．
K_z(t,\tau)=z\left(t+\frac{\tau}{2} \right)z^*\left(t-\frac{\tau}{2} \right)
まずはこれをmatlabで書いてみる．
τの範囲
$\tau$の範囲はウィグナーヴィレ分布の定義では$-\infty$から$+\infty$までだが，有限長の信号だと範囲も有限になる．
分かりやすく，1秒間の信号$z(t),0 \le t \le 1$を考えてみる．
$t=0$の時，
\begin{align}
K_z(0,\tau)&=z\left(0+\frac{\tau}{2} \right)z^*\left(0-\frac{\tau}{2} \right)\\
&=z\left(\frac{\tau}{2} \right)z^*\left(-\frac{\tau}{2} \right)
\end{align}
\left( 0 \le\frac{\tau}{2} \le 1 \ \ \rm{and} \ \ 0 \le -\frac{\tau}{2} \le 1 \right)
となり，$\tau$の取りうる範囲は$\tau=0$のみである．
$\tau>0$の値を取ると$z^*\left( \cdot \right)$の中身が負となって範囲外になってしまうからである．
$t=0.5$の時，
\begin{align}
K_z(0.5,\tau)&=z\left(0.5+\frac{\tau}{2} \right)z^*\left(0.5-\frac{\tau}{2} \right),\\
\end{align}
\left( 0 \le 0.5+\frac{\tau}{2} \le 1 \ \ \rm{and} \ \ 0 \le 0.5-\frac{\tau}{2} \le 1 \right)
となり，$\tau$の取りうる範囲は$-1 \le \tau \le 1$となる．
図にすると以下のような感じとなり，信号の真ん中の部分ほど$\tau$が動ける範囲がでかい．
さらに図にすると↓こんな感じで，値が存在する範囲はひし形となる．
この辺を念頭に置きつつmatlabで書いてみる．
function K_z = mySignalKernel(z)
N = length(z);
K_z = zeros(N);
tau2_center = ceil(N/2);
for t_ind = 1:N
tau2 = min(t_ind - 1, N - t_ind);
for tau2_ind = -tau2:tau2
K_z(t_ind, tau2_ind + tau2_center) ...
= z(t_ind + tau2_ind) * conj(z(t_ind - tau2_ind));
end
end
end
配列のインデックスがぜんぜんわかってないが，雰囲気でこねくり回したら出来た（偶然に基づくプログラミング）．
まぁいいでしょ．ヨシ！
メモリ使用量に注意
ちなみに，さらっと書いたが，
N = length(z);
K_z = zeros(N);
とある通り，信号長さ×信号長さのクソデカ配列が必要となる．
$F_s=44100\rm{[Hz]}$で3分間のWAVファイルなどを取り扱おうとすると…
(8 [Byte] * 44100 [Hz] * 3 [min] * 60 [sec])^2=4.69[GB]？
のメモリが必要なようで（計算が合っているかは知らない），がんばってメモリを増築するか，信号を分割するなど工夫して処理する必要が出てくる．
SignalKernelの可視化
サンプル信号として周波数 $f = 3 \rm{[Hz]}$の$\sin$波を入れて結果を見てみる．
clc
clear all
close all
Fs = 1000; % [Hz]
T = 1; % 信号長さ[s]
L = Fs * T;
t = ((1:L)-1)/Fs;
f1 = 3; % [Hz]
y = sin(2*pi*f1*t);% 周波数一定の信号
z = hilbert(y); % 解析信号にする
% signal kernelを求める
K_z = mySignalKernel(z);
% signal kernelを描く
figure
t = ((1:L)-1)/Fs;
tau = ((1:L)-1)/Fs*2 - L/2/Fs*2;
mesh(t, tau, real(K_z'));
xlabel('t');
ylabel('\tau');
xlim([t(1),t(end)]);
ylim([tau(1),tau(end)]);
view([0,90]);
colormap jet
colorbar
結果はこんな感じ．ひし形に値が入っていることが分かる．
斜めから見た図．
この視点から見ると，確かに周波数 $f = 3 \rm{[Hz]}$ になっていそうだ．
ウィグナーヴィレ分布の導出
\begin{align}
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\\
&=\int ^ \infty _{-\infty} K_z (t, \tau) \mathrm{e}^{-j2 \pi f \tau} d\tau\\
\end{align}
ということで，あとはSignal Kernel $K_z$を$\tau \to f$でフーリエ変換すればおけまる水産なはずである．
function myWVD(K_z, Fs)
L = size(K_z, 1);
% フーリエ変換する
Y = fft(K_z, [], 2);
% この辺はmatlab公式のFFTのサンプルを丸コピペである．
% 心を込めて写経している．つまり何も分かっていない．
P2 = abs(Y/L);
P1 = P2(:,1:L/2+1);
% matlab公式サンプルだと2倍しているが…
%P1(:,2:end-1) = 2*P1(:,2:end-1);
% 今回読み込んだ信号は解析信号なので，
% 片側スペクトルを求める際に2倍する必要はない・・・はず
P1(:,2:end-1) = P1(:,2:end-1);
% グラフを描く
t = ((1:L)-1)/Fs;
f = Fs/2*(0:(L/2))/L;% サンプリング周波数は半分になる
mesh(t,f, P1');
xlabel('t')
ylabel('freq')
view([0,90]);
colormap jet
colorbar
end
ウィグナーヴィレ分布結果
サンプル信号として周波数 $f = 200 \rm{[Hz]}$の$\sin$波を入れて結果を見てみる．
clc
clear all
close all
Fs = 1000; % [Hz]
T = 1; % 信号長さ[s]
L = Fs * T;
t = ((1:L)-1)/Fs;
f1 = 200; % [Hz]
y = sin(2*pi*f1*t);
z = hilbert(y); % 解析信号にする
% signal kernelを求める
K_z = mySignalKernel(z);
% ウィグナーヴィレ分布を求める
myWVD(K_z, Fs);
確かに，$f = 200 \rm{[Hz]}$のところにスペクトルが出ていて，良さそうな感じがする．
ちなみに，サンプリング周波数$F_s=1000\rm{[Hz]}$ならば周波数範囲は$F_s/2=500\rm{[Hz]}$のはずなのに，何故$250\rm{[Hz]}$までしか表示されないかと言うと，Signal Kernelを求める際に，$z(t+\tau/2)z^*(t-\tau/2)$を計算したと思うけど，ここの$\tau/2$に入る値が，$\tau=0,2,4,6,\ldots$と1個飛ばしになっていたからである（サンプリング周波数が半分になる）．じゃないと配列のインデックスに小数点の数値が入っちゃうからね．
チャープ信号のウィグナーヴィレ分布
以下のLFM信号を考える．
f_i(t)=100t+10
% チャープ信号
f1 = 100*t + 10;
% 瞬時周波数のグラフを描く．
figure
plot(t,f1);
xlabel('t');
ylabel('freq');
ylim([0,Fs/4]);
fi = 50*t.^2 + 10*t;
y = sin(2*pi*fi);
% 以降は同様なので省略．
これは求めたい瞬時周波数のグラフ．
ウィグナーヴィレ分布．これもしっかり求められているように見える．
ちなみにSignal Kernelはこんな感じ．横縞模様が時間経過でだんだん細かくなっていっている．
複数信号のウィグナーヴィレ分布
以下の2つのLFM信号を考える．
f_i(t)=100t+10\\
g_i(t)=-30t+150
% チャープ信号
f1 = 200*t + 10;
f2 = -30*t + 150;
figure
plot(t,f1);
hold on
plot(t,f2);
xlabel('t');
ylabel('freq');
ylim([0,Fs/4]);
fi1 = 100*t.^2 + 10*t;
fi2 = - 15*t.^2 + 150*t;
y = sin(2*pi*fi1) + sin(2*pi*fi2);
うむ，こっちもしっかりウィグナーヴィレ分布出来てるんじゃない？そう思うわね．
縞模様が出来ているが，これはクロス項といって，信号と信号の間に発生するスペクトル成分で，ウィグナーヴィレ分布の弱点ね．
ちなみにSignal Kernelはこんな感じ．キモいわね．
これは$(t, \tau)$平面でなく，$(\nu, \tau)$平面で見てみると良さそう～．
今後の課題わね．
matlab標準のwvdの結果
Signal Processing ToolBoxにウィグナーヴィレ分布を求める関数wvdが入っているので，これも表示させてみる．
wvd(z, Fs);
colormap jet
↓wvdの結果．
↓手実装の結果．大体同じ感じになったのでいいんじゃないですかね．
振幅が全く違ってるのがアレなんだよなぁ．振幅$1$の信号を入れてるので，手実装のほうがそれっぽい値になっている感じはするけど…よく分からないですね，雰囲気で信号処理をやっているので…．
まとめ
今回はウィグナーヴィレ分布を手実装してみた．
フーリエ変換をあんまり理解していないのでしっかり出来ているか不安だが，結果自体はそこそこになったので良しとする．
ウィグナーヴィレ分布はSignal Kernel $K_z(t, \tau)$を$\tau \to f$でフーリエ変換すると求められるが，$K_z(t, \tau)$を$\nu \leftarrow t$で逆フーリエ変換すると曖昧度関数$A(\nu, \tau)$が求められる．
クロス項はこの曖昧度関数で見ると分かるようになるらしいので，次回はこの辺をmatlabの実装で結果を確認しつつ数学的にも理解を頑張っていくことにする．
最終的にはウィグナーヴィレ分布を一般化したコーエンクラスをやっていく．