MUSIC 演算法簡介
MUSIC 是一個用於頻率估計以及無線電波方向估計(如 DOA，direction of arrival) 的演算法，全名為 multiple signal classification。MUSIC 在 1980 年代被提出，一開始用於感測器陣列，如雷達、麥克風陣列等。
因為實驗室學長的研究有用到，就來認識一下這個演算法。也順便模擬一下，視覺化呈現 MUSIC。
MUSIC
假設我們知道訊號向量 $\bf x$ 由 $p$ 個 complex exponentials 組成，我們要怎麼估計這 $p$ 個頻率？我們建立一個線性模型：
$$
{\bf x} = A {\bf s} + {\bf n}
$$
其中 $A = [{\bf a_0}, \cdots, {\bf a_{p - 1}}]$，${\bf a_i} = [1, e^{j\omega_i}, e^{j2\omega_i}, \cdots, e^{j(N - 1)\omega_i}]^T$，$\bf s$ 為 complex amplitude。演算法的一個前提為 $p < N$。$\bf x$ 的自相關矩陣為
$$
R_x = A R_s A^H + \sigma^2 I
$$
$R_x = R_x^H$，且是正定矩陣。對 $R_x$ 做特徵值分解 (EVD)，可以得到
$$
\begin{align}
R_x = V \Sigma V^H,\ \ &V = [{\bf v_0}, \cdots, {\bf v_{p - 1}}, {\bf v_p}, \cdots, {\bf v_{N - 1}}]\\[10pt]
&\Sigma =
\begin{bmatrix}
\lambda_0 & \mskip -10mu & \mskip -10mu & \mskip -10mu & \mskip -10mu {\Huge 0} & \mskip -10mu\\[-5pt]
& \mskip -10mu \ddots & \mskip -10mu & \mskip -10mu & \mskip -10mu & \mskip -10mu\\[-5pt]
& \mskip -10mu & \mskip -10mu \lambda_{p - 1} & \mskip -10mu & \mskip -10mu & \mskip -10mu\\[-5pt]
& \mskip -10mu & \mskip -10mu & \mskip -10mu \sigma^2 & \mskip -10mu & \mskip -10mu\\[-5pt]
& \mskip -10mu & \mskip -10mu & \mskip -10mu & \mskip -10mu \ddots & \mskip -10mu\\[-5pt]
& \mskip -10mu {\Huge 0} & \mskip -10mu & \mskip -10mu & \mskip -10mu & \mskip -10mu \sigma^2
\end{bmatrix}_{N\times N}
\end{align}
$$
(特徵值從大到小排列。)前 $p$ 個特徵向量生成 $\cal U_s$，後 $N - p$ 個特徵向量則是和 $\cal U_s$ 正交的雜訊空間 $\cal U_n$。這其實就是主成分分析的概念，前 $p$ 個特徵向量即為功率較大的訊號成分，後 $N - p$ 個特徵向量則是雜訊成分。
令 $U_{\text{noise}} = \left[{\bf v_p}, \cdots, {\bf v_{N - 1}}\right]$，即 $\text{col}(U_{\text{noise}}) = {\cal U_n}$。定義 $d({\bf e})$ 為向量 $\bf e$ 在雜訊空間 $\cal U_n$ 的投影量大小，
$$
d({\bf e}) = \lVert U_{\text{noise}}^H {\bf e}\rVert^2 = \sum_{i = p}^{N - 1}{\lvert {\bf v_i}^H {\bf e}\rvert^2}
$$
如果 ${\bf e} \in {\cal U_s}$，則 $d({\bf e}) = 0$。因此 $d({\bf e})$ 小，表示 $\bf e$ 的訊號成分功率較大、雜訊成分功率較小。反過來說，$1/d({\bf e})$ 愈大，$\bf e$ 的訊號成分功率愈大。要判斷兩個向量 $\bf e_1$ 和 $\bf e_2$ 哪一個更接近於訊號、哪一個更接近於雜訊，只要比較兩者的 $1/d({\bf e})$ 即可。還記得我們一開始想要知道的是訊號的頻率成分嗎？令 ${\bf e} = \left[1, e^{j\omega}, e^{j2\omega}, \cdots, e^{j(N - 1)\omega}\right]^T$，比較頻率軸上任兩點的 $1/d({\bf e})$ 大小，不就可以判斷哪一個頻率更接近於訊號？於是我們定義 MUSIC 的 pseudo-spectrum 為
$$
\hat{P}_{\text{MU}}(e^{j\omega}) = \frac{1}{\lVert U_{\text{noise}}^H {\bf e}\rVert^2}
$$
形成連續的頻譜。根據一開始的假設，訊號由 $p$ 個 complex exponential 組成，只要找到頻譜上數值較高的 $p$ 個點，對應到的 $\omega$ 應該就是組成訊號的 $p$ 個頻率。這就是 MUSIC 演算法。
模擬
網路上有人寫了開源的 MUSIC 函式庫，不過我還是自己寫一遍試試看，畢竟演算法本身相當簡單。圖A為模擬結果，矩陣 $A$ 是由三個頻率的 complex exponential 乘上窗函數 (window function) 組成。窗函數是用來模擬訊號的其他(功率較小的)頻率成分。
結語
MUSIC 演算法概念算是簡單，但是總覺得線性模型的假設和真實情況不太一樣，不過既然已經有許多應用，應該是有實用價值的。
附錄
附錄介紹主成分分析 (principal component analysis，PCA)，複習一下線性代數。假設我們有一個 covariance matrix $R_{xy}$，對它做奇異值分解 (SVD) (這邊都假設 mean 為 zero vector)：
$$
{R_{xy}}_{M\times N} = E\left[{\bf x}{\bf y}^H\right] = U_{M\times K} \Sigma_{K\times K} V^H_{K\times N}
$$
$\Sigma$ 為一對角矩陣，$U$ 的 column vectors 為一 $K$ 維向量空間的一組單位正交基底；$V$ 也是如此。要怎麼解讀這個形式呢？
假設我想要知道 ${\bf s} = A{\bf x}$ 和 ${\bf t} = B{\bf y}$ 的相關性，可以計算它們的 covariance matrix，
$$
R_{st} = E\left[{\bf s}{\bf t}^H\right] = E\left[A{\bf x}{\bf y}^H B^H\right] = AR_{xy}B^H
$$
所以只要知道兩向量的線性變換(矩陣 $A$ 和 $B$)，就可以得知變換後兩向量的 covariance matrix。假設對 ${\bf x}$ 作線性變換，以矩陣表示為 ${\bf x’} = U^H {\bf x}$，對 ${\bf y}$ 作線性變換表示為 ${\bf y’} = V^H {\bf y}$，則變換後兩向量的 covariance matrix 為
$$
R_{x’y’} = U^H R_{xy} V = U^H U \Sigma V^H V = \Sigma
$$
也就是說將 $\bf x$ 投影在 $\text{col}(U) = \text{span}\lbrace{\bf u_0}, \cdots, {\bf u_{K - 1}}\rbrace$，得到座標向量 ${\bf x’} = \left[x_0’, x_1’, \cdots, x_{K - 1}’\right]$；將 $\bf y$ 投影在 $\text{col}(V) = \text{span}\lbrace{\bf v_0}, \cdots, {\bf v_{K - 1}}\rbrace$，得到座標向量 ${\bf y’} = \left[y_0’, y_1’, \cdots, y_{K - 1}’\right]$。這兩個座標向量的相關性可以用 $\Sigma$ 來表示，$y_0’$ 只和 $x_0’$ 有關，$y_1’$ 只和 $x_1’$ 有關，以此類推。
原本 $\bf x$ 和 $\bf y$ 的關係很複雜，$R_{xy}$ 是 full matrix，$y_0$ 和 $x_0 \cdots x_{M - 1}$ 都相關；以座標向量來表示 $\bf x$ 和 $\bf y$ (基底則是一組 latent vector，即主成分)，只需要 $K$ 個共變異數就可以描述兩向量的關係，而且這些關係還是一對一的($x_0$ 愈大，$y_0$ 愈大或愈小)，趨近單變數的函數關係。這就是降維 (dimension reduction) 的概念。推薦觀看李宏毅老師的機器學習教學影片[2]，講解非常精彩。
剛才討論的是 $\bf x$ 和 $\bf y$ 的相關性，那麼如果是 $\bf x$ 和自己的相關性呢？因為 $R_x$ 是半正定矩陣，可以作特徵值分解，$R_x = V \Sigma V^H$。將 $\bf x$ 以 $\lbrace {\bf v_0}, \cdots, {\bf v_{K - 1}}\rbrace$ 為基底表示為座標向量 $\bf x’$，各個座標之間是不相關的。
$$
R_{x’} = V^H R_x V = \Sigma
$$
這有什麼好處？原本我們需要 $M\times M$ 個共變異數來描述資料的組成，以主成分作為基底向量的話，只需要 $K$ 個共變異數和 $K$ 個 latent vector 就可以完整描述該資料，大幅減少所需的儲存空間，同時也將資料變化的程度 (variance) 集中在主成分上：最大的特徵值和對應的主成分，意味著資料在這個主成分的變異量最大；最小的特徵值和對應的主成分，意味著資料在這個主成分的變異量小，既然每一筆資料在該成分的變化不大，資訊量較少，就是比較不重要的成分，可以丟棄。(比如說人臉辨識，假如資料裡每張臉的鼻子都長得差不多，那鼻子這個特徵對於辨識不同人臉就沒有幫助。)
上述的 $R_x$ 和 $R_{xy}$ 都是理論值，實際的情況是我們只有好幾筆資料 $X = \left[{\bf x_0}, {\bf x_1}, \cdots, {\bf x_{n - 1}}\right]$，而不知道 $R_x$。這時候就會用 sample covariance matrix $\hat{R_x}$ 來近似 $R_x$ (資料已經先減去 sample mean)，
$$
\hat{R_x} = \frac{1}{n - 1} X X^H
$$
參考資料
[1] Wikipedia, MUSIC (algorithm), en.wikipedia.org.
Available: https://en.wikipedia.org/wiki/MUSIC_(algorithm)
[2] Hung-yi Lee (Youtube), “ML Lecture 13: Unsupervised Learning - Linear Methods”, youtube.com. Available: https://www.youtube.com/watch?v=iwh5o_M4BNU
[3] Alex (HackMD 筆記), “直觀看待SVD及PCA之關係”, hackmd.io.
Available: https://hackmd.io/@alexmav04/rJG6P7EWO