時頻分析的時頻譜分解

截至目前為止發展的幾種時頻分析方法中,科恩系列分布是最為有力的轉換之一，但是其計算複雜度高於短時距傅立葉變換和小波分析，因此在應用上有所受限。 透過將時頻分布的核分解為不同時頻譜的線性組合，每一份時頻譜都是由簡單的短時距傅立葉變換所計算而得，可以協助我們有效降低分析信號時頻特性的計算量。一般來說當我們使用哈爾小波轉換來分解時頻分部的核的時候，我們可以將STFT的計算量降至只需少量的短時距傅立葉變換即可。

基本概念

根據Cunningham和Williams最早提出的想法，首先我們將時頻分布的核分解為不同短時距傅立葉變換的結果

$C(n\,,k)=\sum _{i=1}^{M}\lambda _{i}|STFT_{i}(n,k)|^{2}$

然後將短時距傅立葉變換的結果以時頻譜的窗函數$w_{i}(n)$和信號$x(n)$加以分解而得

$STFT_{i}(n,k)=\sum _{m=0}^{N-1}x(m-n)w_{i}(m)e^{-j2\pi \,mk/N}$

則我們能順利將核分解為不同時頻譜的窗函數$w_{i}(t)$，但是儘管不同的時頻譜之間彼此正交，不同的時頻譜之間並沒有重要性的順序關聯性讓我們能依據逼近精細度的不同而選擇要保留那些項，因此還是必須計算完整的分解以決定所有的時頻譜。

時頻譜分解

一般來說，一個離散雙線性的時頻分布計算由此時頻分布的核$\psi$所決定，並可表達為內積的形式

$X(n,\omega \,,\psi )=\sum _{n_{1}}\sum _{n_{2}}x(n+n_{1})e^{-j\omega \,(n+n_{1})}\psi (-{\frac {n_{1}+n_{2}}{2}},n_{1}-n_{2})[x(n+n_{2})e^{-j\omega \,(n+n_{2})}]^{*}=<,{\tilde {\psi }}{\bar {S}}_{-n}{\bar {M}}_{-\omega }x,{\bar {S}}_{-n}{\bar {M}}_{-\omega }x>$

其中${\bar {S}}_{-n}$和${\bar {M}}_{-\omega }$分別代表時間位移和頻率位移的算子，而${\tilde {\psi }}$則是代表核$\psi$的作用算子,透過將算子${\tilde {\psi }}$以特徵值分解，我們可以將核分解為線性成分，因而時頻分布本身被拆解為時頻譜的和。具體的來說和可以被分解為

$\sum _{k}\sum _{l}\lambda _{k,l}w_{k}(n_{1})w_{l}(n_{2})$

其中$w_{k}$和$w_{l}$是窗函數，而$\lambda _{k,l}$是特徵值，因此我們重新改寫時頻分析的計算為

$X(n,\omega \,,\psi )=\sum _{n_{1}}x(n+n_{1})e^{-j\omega \,(n+n_{1})}\sum _{n_{2}}[x(n+n_{2})e^{-j\omega \,(n+n_{2})}]^{*}\sum _{k}\sum _{l}\lambda _{k,l}w_{k}(n_{1})w_{l}(n_{2})$

則$\sum _{n_{1}}w_{k}(n_{1})x(n+n1)e^{-j\omega \,(n+n_{1})}$可以被看作信號$x(n)$的短時距傅立葉變換$STFT_{k}(n,\omega )$,因而若令$SP_{k,l}(n,\omega )=STFT_{k}(n,\omega )STFT_{l}^{*}(n,\omega )$

$X(n,\omega \,,\psi )=\sum _{k}\sum _{l}\lambda _{k,l}SP_{k,l}(n,\omega )=\sum _{k}\sum _{l}\lambda _{k,l}STFT_{k}(n,\omega )STFT_{l}^{*}(n,\omega )$

則時頻分布拆解為不同的時頻譜的和，係數$\lambda _{k,l}$代表權重，並可被簡單的由特徵值分解計算而得

降低STFT的計算量

$8\times 8$的離散維格納分布的核矩陣為

${\begin{bmatrix}0&0&0&0&0&0&0&1\\0&0&0&0&0&0&1&0\\0&0&0&0&0&1&0&0\\0&0&0&0&1&0&0&0\\0&0&0&1&0&0&0&0\\0&0&1&0&0&0&0&0\\0&1&0&0&0&0&0&0\\1&0&0&0&0&0&0&0\end{bmatrix}}$

選取8點的哈爾小波作為基底，可計算出特徵值的矩陣

$(\lambda _{k,l})={\begin{bmatrix}1&0&0&0&0&0&0&0\\0&-1&0&0&0&0&0&0\\0&0&0&-1&0&0&0&0\\0&0&-1&0&0&0&0&0\\0&0&0&0&0&0&0&-1\\0&0&0&0&0&0&-1&0\\0&0&0&0&0&-1&0&0\\0&0&0&0&-1&0&0&0\end{bmatrix}}$

可見一開始需要計算$8\times 8=64$次的時頻譜,透過特徵值分解只有少數的非零係數，因此可大大減低計算量。

參考資料

• G.S. Cunningham, W.J. Williams, 「Kernel decomposition of time-frequency distributions」, IEEE Trans. Signal Process. 42 (6) (1994) 1425-1442.
• W.J. Williams, T.-H. Sang, J.C. O』Neill, E.J. Zalubas, 「Wavelet windowed time-frequency distribution decompositions」, in: Proc. SPIE: Advanced Signal Processing Algorithms, Architectures, and Implementations VII, vol. 3162, Soc. of Photo-Optical Instrumentation Engineers, 1997, pp. 149-160.
• Boualem Boashash, 「Time frequency signal analysis and processing」, 2016