MUSIC算法相關(guān)原理知識(shí)（物理解讀+數(shù)學(xué)推導(dǎo)+Matlab代碼實(shí)現(xiàn)）

這篇具有很好參考價(jià)值的文章主要介紹了MUSIC算法相關(guān)原理知識(shí)（物理解讀+數(shù)學(xué)推導(dǎo)+Matlab代碼實(shí)現(xiàn)）。希望對(duì)大家有所幫助。如果存在錯(cuò)誤或未考慮完全的地方，請(qǐng)大家不吝賜教，您也可以點(diǎn)擊"舉報(bào)違法"按鈕提交疑問。

部分來自于網(wǎng)絡(luò)教程，如有侵權(quán)請(qǐng)聯(lián)系本人刪除?

教程鏈接：MUSIC算法的直觀解釋：1，MUSIC算法的背景和基礎(chǔ)知識(shí)_嗶哩嗶哩_bilibili

?MUSIC算法的直觀解釋：2，我對(duì)于MUSIC算法的理解_嗶哩嗶哩_bilibili

https://blog.csdn.net/zhangziju/article/details/100730081

?一、MUSIC算法作用

MUSIC（Multiple Signal Classification），多重信號(hào)分類，是一類空間譜估計(jì)算法。其思想是利用接收數(shù)據(jù)的協(xié)方差矩陣(Rx)進(jìn)行特征分解，分離出信號(hào)子空間和噪聲子空間，利用信號(hào)方向向量與噪聲子空間的正交性來構(gòu)成空間掃描譜，進(jìn)行全域搜索譜峰，從而實(shí)現(xiàn)信號(hào)的參數(shù)估計(jì)。

MUSIC算法通常被用來使用麥克風(fēng)陣列進(jìn)行聲源定位。

例如當(dāng)麥克風(fēng)陣列放在一個(gè)房間中，房間中存在一個(gè)聲源。當(dāng)聲源發(fā)聲時(shí)，陣列會(huì)接收到來自目標(biāo)方向的信號(hào)，但是也會(huì)接受到不同方向的反射信號(hào)。MUSIC算法可以剔除掉其余的反射信號(hào)，選出來自目標(biāo)方向的那個(gè)信號(hào)，從而得到目標(biāo)的方向。

聲波是機(jī)械波，通常利用麥克風(fēng)陣列接受后轉(zhuǎn)換為電信號(hào)進(jìn)行處理。當(dāng)信號(hào)是電磁波時(shí)，例如wifi信號(hào)，我們拿天線陣列進(jìn)行接收，此時(shí)仍然可以利用MUSIC算法把不同方向信號(hào)的角度計(jì)算出來。

二、MUSIC算法原理

MUSIC算法適用于來波為平行波，即目標(biāo)與麥克風(fēng)陣列的距離L遠(yuǎn)大于陣元之間間距d。此時(shí)來自目標(biāo)的信號(hào)相對(duì)于每個(gè)陣元的方位角基本可視為相同。具體如下圖所示：

music算法原理,聲學(xué)相關(guān)算法,算法,matlab

圖1

1.時(shí)延、相位差和目標(biāo)方位角關(guān)系

假設(shè)信號(hào)源發(fā)射信號(hào)為

當(dāng)信號(hào)從聲源目標(biāo)S傳播到陣元1時(shí)，信號(hào)傳播了距離，假設(shè)聲速為c，則耗時(shí)

這樣便會(huì)導(dǎo)致陣元1接受信號(hào)的相位和發(fā)射信號(hào)不一致，會(huì)延遲

那么陣元1最終接收信號(hào)為:

?陣元2相對(duì)于陣元1多走了，那么陣元2接受信號(hào)則為：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

其中，代表陣元2接收信號(hào)和陣元1接收信號(hào)之間的相位差

陣元3相對(duì)于陣元1多走了，那么相應(yīng)的相位差為：

則陣元3接收信號(hào)為：

PS：如果還不好理解因?yàn)槎嘧吡艘欢尉嚯x導(dǎo)致的相位差怎么計(jì)算，可以這樣理解(以陣元1和2為例)：

假設(shè)陣元1接收信號(hào)為

因?yàn)樾盘?hào)到達(dá)陣元2多走了一段距離，那么信號(hào)到達(dá)陣元2的時(shí)間總會(huì)相比陣元1延遲（我們通常稱之為時(shí)延,實(shí)際上相位差就是時(shí)延導(dǎo)致的）

那么陣元2接收信號(hào)則為 $music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

很明顯，和的相位差為 $music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

由數(shù)字信號(hào)處理知識(shí)可以知道，同時(shí)在圖1中陣元1和陣元2的時(shí)延

那么就可以得到相位差

2.MUSIC算法核心原理（思路來源）

MUSIC算法的最終目的：計(jì)算

從以上推導(dǎo)可以知道和兩個(gè)陣元接收信號(hào)的相位差緊密相關(guān)。能求到，就能求得。

那么在理想條件，也就是沒有任何的反射和折射，且只有一個(gè)聲源，這時(shí)直接用兩陣元接收信號(hào)相除就能得到相位差，從而得到目標(biāo)方位角。

但是實(shí)際上會(huì)有很多反射折射信號(hào)被麥克風(fēng)陣列所接收，而且聲源可能不止一個(gè)，此時(shí)該怎么辦呢？這就是MUSIC算法需要解決的問題。

好，那么我們假設(shè)一共有兩個(gè)聲源A,B，發(fā)射信號(hào)分別為和（暫時(shí)不考慮反射和折射）

那么在某時(shí)刻t，三個(gè)陣元接收信號(hào)分別為:

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

那么某段時(shí)間內(nèi)，麥克風(fēng)陣列所接收信號(hào)為:

對(duì)應(yīng)可以寫成:

其中我們已知，是需要求得的，可能已知，可能不知（當(dāng)已知且可逆時(shí)，可以直接用逆矩陣求，但是這樣的情況較少）

能不能通過處理，將等式右邊的消除掉？（MUSIC算法的核心）

怎么處理呢？？？

如果能找到三個(gè)復(fù)數(shù)、和，分別對(duì)3個(gè)陣元接收信號(hào)進(jìn)行幅度和相位變換（用1個(gè)復(fù)數(shù)去乘以一個(gè)信號(hào)，則是對(duì)這個(gè)信號(hào)進(jìn)行了幅度和相位變換），且變換后接收信號(hào)之間完全抵消了，即：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

或用矩陣表示：

把、和分別代入上式可以得到：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

即：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

MUSIC算法在此時(shí)進(jìn)行了一個(gè)假設(shè)，即假設(shè)信號(hào)和信號(hào)是不相關(guān)（MUSIC算法的假設(shè)條件1）的（當(dāng)信號(hào)和信號(hào)線性相關(guān)時(shí)，可以找到一個(gè)非零復(fù)數(shù)，使得）

那么此時(shí)上面公式中和的系數(shù)必須都為0，即：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

上式中系數(shù)和可以直接消除掉，那么可以看出，只要找到、和便可求出和

那么現(xiàn)在問題就轉(zhuǎn)換為了，如何找到這一組復(fù)數(shù)、和？？？

要能找到這一組復(fù)數(shù)，必須滿足：陣元個(gè)數(shù) > 聲源信號(hào)的個(gè)數(shù)（MUSIC算法的假設(shè)條件2）

其實(shí)最終就是解：

和

即，那么MUSIC算法是通過求的最大值（譜峰搜索）來找相應(yīng)的解，對(duì)應(yīng)的也就是相應(yīng)的目標(biāo)方位角。

3.MUSIC算法步驟總結(jié)

窄帶遠(yuǎn)場(chǎng)信號(hào)的DOA數(shù)學(xué)模型為：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

其中X為陣列接收到的信號(hào)矩陣，兩個(gè)維度分別代表：陣元個(gè)數(shù)(number of array elements)、采樣點(diǎn)數(shù)(snapshots)；A為陣列方向矩陣，兩個(gè)維度分別代表：陣元個(gè)數(shù)、信號(hào)方向的方向向量；s為信號(hào)源發(fā)射信號(hào)矩陣，兩個(gè)維度分別代表：信號(hào)源個(gè)數(shù)、采樣點(diǎn)數(shù)；N為噪聲矩陣，兩個(gè)維度分別為陣元個(gè)數(shù)、采樣點(diǎn)數(shù)。

那么陣列接收數(shù)據(jù)的協(xié)方差矩陣為：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

由于信號(hào)和噪聲互相獨(dú)立，數(shù)據(jù)協(xié)方差矩陣可以分解為信號(hào)、噪聲相關(guān)的兩部分，其中Rs是信號(hào)的協(xié)方差矩陣，ARsA^H是信號(hào)部分。

對(duì)R進(jìn)行特征分解有：

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

式中，Us是由R的所有特征值中較大的（信號(hào)源個(gè)數(shù)）個(gè)特征向量組成的子空間，稱為信號(hào)子空間；Un是由R的所有特征值中嬌小的（陣元個(gè)數(shù)-信號(hào)源個(gè)數(shù)）個(gè)特征向量組成的子空間，稱為噪聲子空間。

根據(jù)之前我們所推導(dǎo)的MUSIC算法的條件，要求理想情況下信號(hào)子空間和噪聲子空間正交，也就是信號(hào)子空間中的方向向量a(theta)和噪聲子空間正交：

由于噪聲的存在，實(shí)際上a(theta)和Un并不能完全正交。因此實(shí)際上是通過進(jìn)行最小優(yōu)化搜索來實(shí)現(xiàn)的：

和我們上文所說一樣，MUSIC實(shí)際上是通過譜峰搜索來求最優(yōu)解theta：

PS：由于實(shí)際中陣列接受數(shù)據(jù)是有限的，所以通常由協(xié)方差矩陣的最大似然估計(jì)來代替協(xié)方差矩陣：

總結(jié)以上算法原理，MUSIC算法的步驟為：

1.根據(jù)N個(gè)接收信號(hào)矢量得到下面協(xié)方差矩陣的估計(jì)值：

2.對(duì)第1步得到的協(xié)方差矩陣進(jìn)行特征分解

$music算法原理,聲學(xué)相關(guān)算法,算法,matlab$

3.矩陣會(huì)有M個(gè)特征值。將其從大到小進(jìn)行排列：

其中D(D=信號(hào)源個(gè)數(shù))個(gè)較大的特征值對(duì)應(yīng)信號(hào)，將其對(duì)應(yīng)的特征向量看做信號(hào)部分空間。

M-D(M=陣元個(gè)數(shù))個(gè)較小的特征值對(duì)應(yīng)噪聲，將其對(duì)應(yīng)的特征向量看做信號(hào)部分空間，得到噪聲矩陣

4.使不斷變化，計(jì)算譜函數(shù)：

通過尋找譜峰來計(jì)算波達(dá)方向的估計(jì)值。此處的為陣元的方向相應(yīng)向量。

，

4.Matlab代碼實(shí)現(xiàn)文章來源地址http://www.zghlxwxcb.cn/news/detail-787864.html

clear all
close all
clc
%----------------均勻線列陣實(shí)現(xiàn)MUSIC算法------------------%
ang2rad = pi/180;                   % 角度轉(zhuǎn)弧度系數(shù)
N = 10;                             % 陣元個(gè)數(shù)
M = 3;                              % 信源個(gè)數(shù)
theta = [-65,0,45];                 % 來波方向(角度)
snr = 10;                           % 信號(hào)信噪比dB
K = 512;                            % 總采樣點(diǎn)
delta_d = 0.05;                     % 陣元間距
f = 2400;                           % 信號(hào)源頻率
c = 340;                            % 聲速

d = 0:delta_d:(N-1)*delta_d;
A = exp(-1i*2*pi*(f/c)*d.'*sin(theta*ang2rad));   % 接收信號(hào)方向向量
S = randn(M,K);                     % 陣列接收到來自聲源的信號(hào)
X = A*S;                            % 最終接收信號(hào)，是帶有方向向量的信號(hào)
X1 = awgn(X,snr,'measured');        % 在信號(hào)中添加高斯噪聲
Rx = X1*X1'/K;                      % 協(xié)方差矩陣
[Ev,D] = eig(Rx);                   % 特征值分解
% [V,D] = eig(A) 返回特征值的對(duì)角矩陣 D 和矩陣 V
% 其列是對(duì)應(yīng)的右特征向量，使得 AV = VD
EVA = diag(D)';                     % 將特征值提取為1行
[EVA,I] = sort(EVA);                % 對(duì)特征值排序,從小到大。其中I為index：1,2，...，10
EV = fliplr(Ev(:,I));               % 對(duì)應(yīng)特征矢量排序
 En = EV(:,M+1:N);                  % 取特征向量矩陣的第M+1到N列特征向量組成噪聲子空間
 
% 遍歷所有角度，計(jì)算空間譜
for i = 1:361
    angle(i) = (i-181)/2;           % 映射到-90度到90度
    theta_m = angle(i)*ang2rad;
    a = exp(-1i*2*pi*(f/c)*d*sin(theta_m)).';
    p_music(i) = abs(1/(a'*En*En'*a));
end
p_max = max(p_music);
p_music = 10*log10(p_music/p_max);  % 歸一化處理
figure()
plot(angle,p_music,'b-')
grid on
xlabel('入射角/度')
ylabel('空間譜/dB')

到了這里，關(guān)于MUSIC算法相關(guān)原理知識(shí)（物理解讀+數(shù)學(xué)推導(dǎo)+Matlab代碼實(shí)現(xiàn)）的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！