MUSIC算法

**

MUSIC算法仿真分析
一、DOA估计方法
波达方向(Direction of arrival, DOA)是指空间信号的到达方向(各个信号到达阵列参考阵元的方向角，简称波达方向)。它是空间谱估计理论的重要基础，而空间谱估计又是阵列信号处理的研究领域。空间谱估计主要研究空间多传感器阵列所构成的处理系统对感兴趣的空间信号的多种参数进行准确估计的能力 ，主要目的是估计信号的空域参数或信源位置，而波达方向(DOA)正是其中最主要的一个参数之一，对于远场信号，当信号到达不同的阵元时存在一个波程差，波程差导致了接收阵元的到达时间差（time difference of arrival，TDOA），TDOA是一种无源定位技术。

图1基于TDOA 的传感器阵列定向示意图 图2等距线阵
对于图1设K为声源信号，则其单位向量为：

一、MUSIC基本原理
多重信号分类( Multiple Signal Classification，MUSIC)算法的基本原理是把阵列输出数据的协方差矩阵进行特征分解，得到与信号分量相对应的信号子空间和与信号分量正交的噪声子空间，然后利用两个子空间的正交性实现信号的入射方向估计。
设空间有 D 个互不相关的窄带远场信号以方位角 θ1,θ2,…,θD 入射到测向阵列中，入射信号的数目 D 小于阵列的阵元数 M。时间由第 k 个采样时钟表示，则阵列的输出矢量为：
(1)
其中，为阵列接收信号，为远场入射信号，为阵列的方向矩阵，其中为的方向阵列，为测量噪声：
(2)
(3)
(4)
(5)
其中， (6)
(7)
阵列输出矢量的协方差矩阵为：
(8)
其中，为对称矩阵；
(9)
(10)
为信号协方差矩阵，对进行特征值分解：
(11)
在入射信号互不相关的情况下，将的特征值按大小顺序排列，可得到：
(12)
前个特征值与信号有关，其特征值均大于，所以信号子空间由个较大特征值对应的特征矢量构成其余个较小特征值对应的矢量构成噪声子空间,在理想条件下信号子空间和噪声子空间正交，因此有
(13)
实际应用中，根据次快拍得到接收数据，用时间平均估计空间相关，为
(14)
所以MUSIC功率谱函数为：
， (15)
按方位角度进行搜索，取得峰值就是个信号源的波达估计值。
MUSIC算法的计算步骤：

（1）利用阵列所接收的L快拍数据，可以得到数据的协方差矩阵；
（2）对进行特征值分解，由的特征值就可以判断信号源数；
（3）得到信号子空间与噪声子空间；
（4）根据式子进行谱峰搜索；
（5）找出极大值所对应的角度就是入射信号的角度。
二，music算法程序
clear all
close all
clc
%basic parameter

p=3; %信源数
snr=-10; %信噪比
inr=10;
N=1000; %采样快拍数
lamda=1; %信号波长
d=lamda/2; %阵元间距
L=8*lamda; %孔径
m=64; %阵元数
%pre-allocate variables for speed
mag_A=zeros(m,p);sig=zeros(p,N);noise=zeros(m,N);
rx=zeros(m,m);x=zeros(m,N);%定义空矩阵

theta=[0 2 40]; %p个信源的波达方向
a_theta1=(exp(-j2pi*[0:d:(m-1)d]sin(theta(1)pi/180))/lamda)’;
a_theta2=(exp(-j2pi[0:d:(m-1)d]sin(theta(2)pi/180))/lamda)’;
a_theta3=(exp(-j2pi[0:d:(m-1)*d]*sin(theta(3)*pi/180))/lamda)’;
mag_A=[a_theta1,a_theta2,a_theta3]; %m个阵元接收到p个信源的方向引导矩阵

%create gaussian white noise
noise=(randn(m,N)+j*randn(m,N))/sqrt(2);

%create input signal
amp1=sqrt(cov(noise(1,:))*10^(snr/10)); %信号源复包络的幅度值
amp2=sqrt(cov(noise(1,:))*10^(inr/10));
amp3=sqrt(cov(noise(1,:))*10^(inr/10));

f=4000; %频率
f1=500; f2=600;f3=700; %信号频率
t=[0:1/f:(N-1)/f];
sig1=amp1sin(2pif1t);sig2=amp2sin(2pif2t);sig3=amp3sin(2pif3t); %m个阵元接收到p个信源发出的复包络信号
sig=[sig1;sig2;sig3];

%create receive signal
x=mag_Asig+noise; %m阵元的第k次快拍的采样值 k=1,2,…,N
rs=(1/N)(sigsig’);%%
rx=(1/N)(x*x’); %求得自相关函数rx的值

[V,D]=eigs(rx,m-p,‘SM’); %特征值分解
theta0=[-90:0.1:90];
for n=1:length(theta0);
a=(exp(-j2pi*[0:d:(m-1)*d]*sin(theta0(n)*pi/180))/lamda)’;
beam(n)=1/(a’VV’*a); %谱峰搜索
end

p=3; %信源数
snr=-10; %信噪比
inr=10;
N=1000; %采样快拍数
lamda=0.5; %信号波长
d=lamda/2; %阵元间距
L=8*lamda; %孔径
m=64; %阵元数
%pre-allocate variables for speed
mag_A=zeros(m,p);sig=zeros(p,N);noise=zeros(m,N);
rx=zeros(m,m);x=zeros(m,N);%定义空矩阵

theta=[0 2 40]; %p个信源的波达方向
a_theta1=(exp(-j2pi*[0:d:(m-1)d]sin(theta(1)pi/180))/lamda)’;
a_theta2=(exp(-j2pi[0:d:(m-1)d]sin(theta(2)pi/180))/lamda)’;
a_theta3=(exp(-j2pi[0:d:(m-1)*d]*sin(theta(3)*pi/180))/lamda)’;
mag_A=[a_theta1,a_theta2,a_theta3]; %m个阵元接收到p个信源的方向引导矩阵

%create gaussian white noise
noise=(randn(m,N)+j*randn(m,N))/sqrt(2);

%create input signal
amp1=sqrt(cov(noise(1,:))*10^(snr/10)); %信号源复包络的幅度值
amp2=sqrt(cov(noise(1,:))*10^(inr/10));
amp3=sqrt(cov(noise(1,:))*10^(inr/10));

f=4000; %频率
f1=500; f2=600;f3=700; %信号频率
t=[0:1/f:(N-1)/f];
sig1=amp1sin(2pif1t);sig2=amp2sin(2pif2t);sig3=amp3sin(2pif3t); %m个阵元接收到p个信源发出的复包络信号
sig=[sig1;sig2;sig3];

%create receive signal
x=mag_Asig+noise; %m阵元的第k次快拍的采样值 k=1,2,…,N
rs=(1/N)(sigsig’);%%
rx=(1/N)(x*x’); %求得自相关函数rx的值

[V,D]=eigs(rx,m-p,‘SM’); %特征值分解
theta0=[-90:0.1:90];
for n=1:length(theta0);
a=(exp(-j2pi*[0:d:(m-1)*d]*sin(theta0(n)*pi/180))/lamda)’;
beam1(n)=1/(a’VV’a); %谱峰搜索
end
plot(theta0,10log10(abs(beam)/max(beam)),‘r-’);xlabel(‘俯仰角(degree)’);ylabel(‘信噪比（dB）’);grid on;
title(‘经典MUSIC算法DOA估计，(mx，2018042102)’);
hold on

plot(theta0,10*log10(abs(beam1)/max(beam1)),‘r-’);