一、簡介

MUSIC演算法 是一種基於矩陣特徵空間分解的方法。從幾何角度講，訊號處理的觀測空間可以分解為訊號子空間和噪聲子空間，顯然這兩個空間是正交的。訊號子空間由陣列接收到的資料協方差矩陣中與訊號對應的特徵向量組成，噪聲子空間則由協方差矩陣中所有最小特徵值(噪聲方差)對應的特徵向量組成。
1 演算法原理
MUSIC演算法是空間譜估計測向理論的重要基石。如下:

(1) 不管測向天線陣列形狀如何，也不管入射來波入射角的維數如何，假定陣列由M個陣元組成，則陣列輸出模型的矩陣形式都可以表示為:Y(t)=AX(t)+N(t)

其中，Y是觀測到的陣列輸出資料復向量;X是未知的空間訊號復向量;N是陣列輸出向量中的加性噪聲;A是陣列的方向矩陣;此處，A矩陣表示式由圖冊表示。

MUSIC演算法的處理任務就是設法估計出入射到陣列的空間訊號的個數D以及空間訊號源的強度及其來波方向。

(2) 在實際處理中，Y得到的資料是有限時間段內的有限次數的樣本(也稱快拍或快攝)，在這段時間內，假定來波方向不發生變化，且噪聲為與訊號不相關的白噪聲，則定義陣列輸出訊號的二階矩:Ry。

(3) MUSIC演算法的核心就是對Ry進行特徵值分解，利用特徵向量構建兩個正交的子空間，即訊號子空間和噪聲子空間。對Ry進行特徵分解，即是使得圖冊中的公式成立。

(4) U是非負定的厄米特矩陣，所以特徵分解得到的特徵值均為非負實數，有D個大的特徵值和M-D個小的特徵值，大特徵值對應的特徵向量組成的空間Us為訊號子空間，小特徵值對應的特徵向量組成的空間Un為噪聲子空間。

(5) 將噪聲特徵向量作為列向量，組成噪聲特徵矩陣 ，並張成M-D維的噪聲子空間Un，噪聲子空間與訊號子空間正交。而Us的列空間向量恰與訊號子空間重合，所以Us的列向量與噪聲子空間也是正交的，由此，可以構造空間譜函式。

(6) 在空間譜域求取譜函式最大值，其譜峰對應的角度即是來波方向角的估計值。

2 理論發展及應用
MUSIC(Multiple Signal Classification多訊號分類)演算法是1979年由美國人R.O.Schmidt提出的，它標誌著空間譜估計測向進入了繁榮發展的階段。它將"向量空間"的概念引入了空間譜估計領域，經過三十年的發展，可以說其理論已經比較成熟。

自80年代以來，人們對基於特徵分解的超解析度空間譜估計演算法進行了廣泛深入的研究，並提出了一系列高效的處理方法，其中最經典的是多訊號分類(MUSIC)演算法，這種演算法要經過一維搜尋才能求出信源的來向，而相對最大似然(ML)和加權子空間擬合(WSF)等多維搜尋演算法的運算量已經減少了很多。以MUSIC為代表的演算法存在一個缺點，即對相干訊號處理的不理想。在針對相干訊號源的一系列處理方案中，比較經典的是空間平滑技術，如空間平滑(SS)和修正的空間平滑(MSS)演算法。然而，空間平滑技術是以損失陣列有效孔徑為代價的，而且只適用於等距均勻線陣(ULA)。

事實上空間譜估計演算法都是在已知訊號源數目下計算的，而在實際應用中這是不可能的，只能根據觀測資料對源數目進行估計。R.O.Schmidt在他的經典之作中提出了依據陣列協方差矩陣特徵值的分佈來估計訊號源的方法。這種方法在理論上是完美的，至少對獨立源和部分相關源是正確的，但實際上由於資料長度有限，很大程度上只能依靠主觀判斷來確定源數。

二、原始碼

clc;
clear all;
close all;
set(0,'defaultaxesfontsize',9);            %設定字型大小
load s.mat;
wnd=256;
inc=128;

[Angle1]=Spectrum_Method('capon',s1,s2,wnd,inc,45);
subplot(311)
plot(Angle1,'*'),axis tight
title('Capon')
xlabel('幀數')
ylabel('角度/度')

[Angle2]=Spectrum_Method('music',s1,s2,wnd,inc,45);
subplot(312)
plot(Angle2,'*'),axis tight
title('Music')
xlabel('幀數')
ylabel('角度/度')
function frameout=enframe(x,win,inc)

nx=length(x(:));            % 取資料長度
nwin=length(win);           % 取窗長
if (nwin == 1)              % 判斷窗長是否為1，若為1，即表示沒有設窗函式
   len = win;               % 是，幀長=win
else
   len = nwin;              % 否，幀長=窗長
end
if (nargin < 3)             % 如果只有兩個引數，設幀inc=幀長
   inc = len;
end
nf = fix((nx-len+inc)/inc); % 計算幀數
frameout=zeros(nf,len);            % 初始化
indf= inc*(0:(nf-1)).';     % 設定每幀在x中的位移量位置
inds = (1:len);             % 每幀資料對應1:len
frameout(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:));   % 對資料分幀

三、執行結果

四、備註

版本：2014a
完整程式碼或代寫加1564658423