重新分布法

重新分布法（Reassignment Method）是一種銳化訊號之時頻表示的方法，藉由將資料映射至較靠近原始訊號之真實支撐區域的時頻座標來實現。此方法曾被不同學者獨立提出，並有重映射（remapping）、時頻重新分布（time-frequency reassignment）以及修正滑動視窗法（modified moving-window method）等別稱。以時頻譜或短時距傅立葉變換而言，重新分布法可藉由估算局部的瞬時頻率以及群延遲，使模糊的時頻資料點重新定位並清晰化。當訊號可藉由分析視窗進行時域和頻域的分離時，這項時頻座標的重新映射是相當精準的。

## 概述

<image>

許多信號的能量分布都具有隨著時間和頻率變化的特性。舉例而言，任何有起始或者結束的音訊訊號之能量分布都將隨著時間改變，而大多數聲音的能量不論在時間上還是頻率上都有相當大的變化；時頻分析便經常用於尋找此類訊號的特徵，這種分析法將原為一維的時域訊號轉換為二維的時域–頻域函數。而經分析轉換後的結果便稱為時頻表示，它可以用於描述一個訊號隨時間的頻譜能量分布變化，就像樂譜描述樂曲中高低音隨時間的變化一樣。

在音訊訊號分析中，時頻譜是最常用的時頻表示法；這可能是因為該方法長久以來為學界所熟知，且能夠避免產生不易解讀的「交叉項」（cross-terms）。不過，計算時頻譜所須使用的視窗函數使得分析結果必須在時域或頻域的解析度中做一取捨，也就是時頻譜在時域及頻域的解析度必定至少有其中一者不理想。這項時頻表示的重新分布法便是透過先把資料映射至較接近原始訊號之真實支撐區域的座標，以將時頻譜或其它時頻表示下較模糊的資料點重新定位並聚焦。

## 時頻譜

時頻譜為最為人所熟知的時頻表示法之一，其定義為短時距傅立葉變換值的平方：

$$SP_x(t,f)=|X(t,f)|^2=\left|\int_{-\infty}^{\infty} w(t-\tau)x(\tau)e^{-j2\pi f \tau}d\tau\right|^2.$$

雖然短時距的相位頻譜確實含有與訊號相關的重要時態資訊，但這些資訊通常難以解讀，故一般而言在短時距頻譜分析中只考慮其短時距強度頻譜。以時頻表示法而言，時頻譜具有相對較差的解析度。它在時域與頻域的解析度取決於分析時所選擇的視窗函數，且其中一者愈清晰將導致另一者愈模糊。

相對於時頻譜，維格納–魏爾分布則是另一種解析度較高的時頻表示法；本質上，維格納–魏爾分布的做法與短時距傅立葉變換類似，惟其選擇的視窗函數與訊號之間有完好的對應。維格納–魏爾分布在時域和頻域中皆可得到高解析度，但它同時也具有高度的非線性和非局部性等特質，因此對雜訊相當敏感且容易產生覆蓋訊號本身的交叉項，使得含有多成分的訊號分析後的資訊難以解讀。

利用雙線性轉換得到時頻表示的科恩系列分布函數則可視為「平滑化」的維格納–魏爾分布。它藉由一平滑核降低分布對於雜訊的敏感度，並抑制交叉項的產生，但此舉將使時域和頻域的分布模糊化；這不僅會降低頻譜的解析度，同時也可能使原始維格納準機率分布中實際上沒有能量的區域得到非零的能量值。

時頻譜便屬於科恩系列分布之一；它也可視為一種平滑化的維格納–魏爾分布，惟其使用之平滑核即相等於其分析視窗函數的維格納–魏爾分布。重新分布法便是先將維格納–魏爾分布平滑化，再將它於原始訊號成分之真實支撐區域附近的座標點重新聚焦，進而達到銳化的目的。實驗結果顯示，此方法可以減少任一科恩系列分布在時域與頻域的模糊現象。對時頻譜而言，重新分布法便是利用訊號的短時距相位頻譜來修正其時頻座標點，再將其重新映射回較靠近所分析訊號之真實支撐區域的位置。

## 方法

重新分布法首先由小寺邦彥（Kunihiko Kodera）、羅傑·干德林（Roger Gendrin）和克勞地·德·凡德瑞（Claude de Villedary）三人提出，並命名為「修正滑動視窗法」（Modified Moving-Window Method）。他們的方法為每個資料點指定一組新的時頻座標，使其更佳反映原始訊號的能量分布，進而改善了標準滑動視窗法（等同於時頻譜）的時頻解析度。

標準的滑動視窗法乃是將一時域訊號$x(t)$分解為一組基底訊號$h_{\omega}(t)$的加權總和，其中代表每項基底的係數則為$\epsilon( t, \omega )$。基底訊號的定義如下：

$$h_{\omega}(t) = h(t) e^{j \omega t}$$

其中，$h(t)$為一實值之低通核函數，如短時距傅立葉變換中的視窗函數。此般分解下的各項係數則可定義如下：

$$\begin{align} \epsilon( t, \omega ) &= \int x(\tau) h( t - \tau ) e^{ -j \omega \left[ \tau - t \right]} d\tau \\ <pre> &= e^{ j \omega t} \int x(\tau) h( t - \tau ) e^{ -j \omega \tau } d\tau \\ &= e^{ j \omega t} X(t, \omega) \\ &= X_{t}(\omega) \\ &= M_{t}(\omega) e^{j \phi_{\tau}(\omega)}. </pre> \end{align}$$

其中，$X_{t}(\omega)$為訊號$x(t)$在時間位移為$t$、使用$h(t)$做為視窗函數下所得到的傅立葉變換；而$M_{t}(\omega)$與$\phi_{\tau}(\omega)$則分別代表$X_{t}(\omega)$的強度以及相位值。

$x(t)$可透過滑動視窗的係數以下式重建：

$$\begin{align} x(t) & = \iint X_{\tau}(\omega) h^{*}_{\omega}(\tau - t) d\omega d\tau \\ <pre> & = \iint X_{\tau}(\omega) h( \tau - t ) e^{ -j \omega \left[ \tau - t \right]} d\omega d\tau \\ &= \iint M_{\tau}(\omega) e^{j \phi_{\tau}(\omega)} h( \tau - t ) e^{ -j \omega \left[ \tau - t \right]} d\omega d\tau \\ &= \iint M_{\tau}(\omega) h( \tau - t ) e^{ j \left[ \phi_{\tau}(\omega) - \omega \tau+ \omega t \right] } d\omega d\tau. </pre> \end{align}$$

當一個訊號的強度頻譜為$M(t,\omega)$，且其時間上較相位上的變化來得慢時，上述重建式中積分值最主要的貢獻來源即滿足下列穩相條件（phase stationarity condition）的$t,\omega$及其鄰近區域：

$$\begin{align} \frac{\partial}{\partial \omega} \left[ \phi_{\tau}(\omega) - \omega \tau + \omega t\right] & = 0, \\ \frac{\partial}{\partial \tau} \left[ \phi_{\tau}(\omega) - \omega \tau + \omega t \right] & = 0. \end{align}$$

亦即滿足下列定義的$\hat{t}, \hat{\omega}$及其鄰近區域：

$$\begin{align} \hat{t}(\tau, \omega) & = \tau - \frac{\partial \phi_{\tau}(\omega)}{\partial \omega} = -\frac{\partial \phi(\tau, \omega)}{\partial \omega}, \\ \hat{\omega}(\tau, \omega) & = \frac{\partial \phi_{\tau}(\omega)}{\partial \tau} = \omega + \frac{\partial \phi(\tau, \omega)}{\partial \tau}. \end{align}$$

此現象在光學等領域中稱為穩相近似：對週期或準週期函數而言，其傅立葉相位頻譜在其振盪頻率附近相對於時間的變化較慢，反之在外圍區域則較快。以時域上能量較為集中的脈衝訊號為例，其相位頻譜在脈衝時間點附近的變化較慢，而較遠的區域則變化較快。

在訊號重建時，合成波形之正值與負值於相位快速變化的頻率區域會因破壞性干涉而互相抵消。只有相位變化緩慢（即穩相）的區域會對重建的訊號有顯著的影響，故貢獻極大值（訊號的重心）即發生在相位隨時間與頻率變化最緩慢的地方。

由此知，可以直接從（一般在建構時頻譜時被忽略的）短時距傅立葉轉換之相位值求得時頻座標，且其值即等於訊號的局部群延遲$\hat{t}_g(t,\omega)$和瞬時頻率$\hat{\omega}_i(t,\omega)$。這些時頻值代表的是經由視窗函數過濾後的訊號，因而屬於局部（local）的量且無法表示整個訊號的全域（global）特徵。

修正滑動視窗法（亦即重新分布法）便是將計算$\epsilon(t,\omega)$的時頻點從原始的$t,\omega$重新移至局部極大貢獻點$\hat{t}(t,\omega), \hat{\omega}(t,\omega)$。由於本質上這個局部極大點和質量分布中心點有相似之處，故也常類比成「重心」。該類比也有助於做出下述的理解：一個訊號的頻譜能量分布對重心的貢獻只在「有能量」的區域才具有意義，因此對時頻譜上量值為零的座標點使用重新分布法是沒有意義的。

## 快速時頻重新分布法

在數位訊號處理中，最常見的做法是先在時域與頻域分別對訊號進行取樣，並使用離散傅立葉變換計算時域訊號$x(n)$的頻域變換$X(k)$。因離散的時頻資料點無法直接計算偏微分，小寺邦彥等人所提出的重新分布法並不能直接套用於離散形式的短時距傅立葉變換；甚至有人指出，這項難處就是使得重新分布法無法得到更廣泛應用的主因。

然而，偏微分的值可以利用有限差分法來估計。舉例而言，在相位頻譜上兩個鄰近的時間點（$t+\Delta t/ 2$、$t-\Delta t/ 2$）取值、並將兩值相減後除以時間差$\Delta t$，便可得到該處頻譜於時域的偏微分估計值（頻域亦同理）：

$$\begin{align} \frac{\partial \phi(t, \omega)}{\partial t} & \approx \frac{1}{\Delta t} \left[ \phi \left (t + \frac{\Delta t}{2}, \omega \right ) - \phi \left (t - \frac{\Delta t}{2}, \omega \right ) \right] \\ \frac{\partial \phi(t, \omega)}{\partial \omega} & \approx \frac{1}{\Delta \omega}\left[ \phi \left (t, \omega+ \frac{\Delta \omega}{2} \right ) - \phi \left (t, \omega-\frac{\Delta \omega}{2} \right ) \right] \end{align}$$

在頻譜上，某些區域的相位變化主要源自於附近某單一弦波成分的旋轉；這些區域內的相位值為一線性函數，因此對於足夠微小的$\Delta t$及$\Delta \omega$而言，上式的有限差分可為相位值的偏微分做一良好的近似。

道格拉斯·尼爾森（Douglas Nelson）也曾獨立於小寺邦彥等人，開發了利用短時距相位頻譜的偏微分來改善時頻譜之解析度的類似方法。稍加推導便可容易地證明他提出的「跨頻譜平面」（cross-spectral surface）所計算之偏微分近似值和前述的有限差分法無異。

弗朗索瓦·俄歇（François Auger）及派翠克·弗朗德蘭（Patrick Flandrin）則指出，小寺邦彥等人所提出的時頻譜重新分布法可透過下述的一般化，延用於任一科恩系列的時頻表示法：

$$\begin{align} \hat{t} (t,\omega) &= t - \frac{\iint \tau \cdot W_{x}(t-\tau,\omega -\nu) \cdot \Phi(\tau,\nu) d\tau d\nu} {\iint W_{x} \left (t-\tau,\omega -\nu \right ) \cdot \Phi (\tau,\nu) d\tau d\nu } \\ \hat{\omega} (t,\omega) & = \omega - \frac{\iint \nu \cdot W_{x}(t-\tau,\omega -\nu) \cdot \Phi(\tau,\nu) d\tau d\nu} {\iint W_{x}(t-\tau,\omega -\nu) \cdot \Phi(\tau,\nu) d\tau d\nu} \end{align}$$

其中，$W_{x}(t,\omega)$為$x(t)$的維格納–魏爾分布，而$\Phi(t,\omega)$則是定義該分布的核函數。他們更提出了一項不須直接計算相位之偏微分，便可快速且準確地得到重新分布後之時頻譜的方法。以時頻譜而言，其重新分布法的計算如下：

$$\begin{align} <pre> \hat{t} (t,\omega) & = t - \Re \left \{ \frac{ X_{\mathcal{T}h}(t,\omega) \cdot X^*(t,\omega) }{ | X(t,\omega) |^2 } \right \} \\ \hat{\omega}(t,\omega) & = \omega + \Im \left \{ \frac{ X_{\mathcal{D}h}(t,\omega) \cdot X^*(t,\omega) }{ | X(t,\omega) |^2 } \right \} </pre> \end{align}$$

其中，$X(t,\omega)$是訊號以視窗函數$h(t)$所得到的短時距傅立葉變換、$X_{\mathcal{T}h}(t,\omega)$是訊號以時間加權之視窗函數$h_{\mathcal{T}}(t) = t \cdot h(t)$所得到的短時距傅立葉變換、$X_{\mathcal{D}h}(t,\omega)$則是訊號以對時間取導數之視窗函數$h_{\mathcal{D}}(t) = \tfrac{d}{dt}h(t)$所得到的短時距傅立葉變換；故知，利用輔助視窗函數$h_{\mathcal{T}}(t)$以及$h_{\mathcal{D}}(t)$，可得到含原始在內共三種傅立葉變換的結果。接著，將這三種變換分別在$t,\omega$處取值並做代數上的組合，便可在任一時頻座標點$t,\omega$計算重新分布法的結果。由於上述的方法只是在短時距頻譜上的某單一時頻點取值，而不須做任何導數的計算，因此可視作離散短時距傅立葉變換之重新分布的快速演演算法。

上述方法有一限制，即$| X(t,\omega) |^2$之值必須為非零；不過，這項限制一般而言皆會被滿足，因為只在頻譜有能量的區域才有值可以「重新分布」，而方法本身在頻譜值為零處也沒有意義。

## 可分離性

<image>

<image>

短時距傅立葉變換經常用於估算一個多成分訊號（例如準諧波的樂器音）中個別單一成分的振幅和相位，而時頻重新分布則能將變換後的時頻點集中至能量分布的局部重心，使頻譜更加清晰。

對單一成分的訊號而言，其瞬時頻率可由該訊號經任意短時距傅立葉變換後得到的相位之偏微分估算。若訊號可拆解為多個成分：

$$x(t) = \sum_{n} A_{n}(t) e^{j \theta_{n}(t)}$$

且定義各成分的瞬時頻率為其相位對時間的偏微分，即：

$$\omega_{n}(t) = \frac{d \theta_{n}(t)}{d t}$$

則可以藉由設計適當的帶通濾波器使其通帶只包含某一特定訊號成分，再從該濾波器經變換後之頻譜的相位值，計算該訊號成分的瞬時頻率。

尼爾森稱這項頻域中的條件為「可分離性」（separability），且所有透過此方法分析的訊號皆須滿足該條件。若該條件並未滿足，則因不能透過短時距傅立葉變換估算個別成分的參數，而無法得到欲求的多成分訊號拆解；在這樣的情況下，必須採用其它的分析視窗函數，以滿足可分離性的條件。

若某特定的短時距頻譜分析視窗可用於拆解訊號中不同頻率的成分，則每個濾波器經短時距傅立葉變換後的脈衝響應便是（至多一項）單一主要訊號成分（即能量最大）經濾波後的結果，故$X(t,\omega_0)$之相位值對時間的偏微分值，即相等於頻率$\omega_0$下主要成分之相位值對時間的偏微分值。因此，若某瞬時頻率為$\omega_{n}(t)$的訊號成分$x_n(t)$ 在$\omega_0$附近為主要成分時，則該成分的瞬時頻率可透過原始訊號經短時距傅立葉變換後位於$\omega_0$的相位值求得，即：

$$\begin{align} \omega_{n}(t) &= \frac{\partial}{\partial t} \arg\{ x_{n}(t) \} \\ <pre> &= \frac{\partial }{\partial t} \arg\{ X(t,\omega_{0}) \} </pre> \end{align}$$

頻域的可分離性要求每個帶通濾波器經短時距傅立葉變換後，最多只可通過一項訊號成分，而在時域同樣也有類似的限制：訊號中兩個不同時間的成分（即事件）彼此間分隔須夠大，使得任一時域視窗下的分段，都不得包含超過一個事件。此即時域下的可分離性條件，且等義於任兩事件間的間隔須大於短時距傅立葉變換濾波器之脈衝響應長度，也就是$h(t)$中非零取樣點所涵蓋的範圍。

一般而言，任一多成分訊號有無限多種相同合理的拆解方式；不過，可分離性必須是在所欲求得的拆解下進行考慮。舉例而言，在分析一個語音訊號時，若使用時間範圍較聲門脈衝間隔長的分析視窗，則足以將不同諧波分離；但個別的聲門脈衝則會因單一視窗包含了多個脈衝而被模糊，也就是所選的分析視窗使得訊號在時域中不具可分離性。反之，若使用時間範圍遠小於聲門脈衝間隔的分析視窗，則能夠成功分離每個脈衝，因為每個視窗中最多只包含一個脈衝；但不同諧波便會彼此模糊交疊，因為視窗在頻域中的範圍比單一諧波還要寬，也就是所選的分析視窗使得訊號在頻域中不具可分離性。

## 外部連結

- TFTB — Time-Frequency ToolBox
- SPEAR - Sinusoidal Partial Editing Analysis and Resynthesis
- Loris - Open-source software for sound modeling and morphing
- SRA - A web-based research tool for spectral and roughness analysis of sound signals
- Sparse time-frequency representations - PNAS