untitled.m

利用music方法来进行doa估计

%%信号为矩形单频信号，噪声为高斯白噪声，阵列为8元均匀线阵列
%%用music算法完成两不相关信号源的方位估计；信噪比分别为：-5dB,0dB,5dB,10dB.

%文件头
clear all;
close all
clc;
tic;
%第一步： 设置参数
 M=8;                    %阵元数
 N=2;                    %信号源数目
 N_s=1024;               %信号长度
 l=0.15;                 %波长，                                              ！！??此处的参数设置有些不符合实际？
 d=0.5*l;                %阵元间隔，为半波长
 SNR=10;                 %信噪比
 w=2*pi/16;              %信号频率
 j=sqrt(-1);
 xx=-pi/6;  yy=pi/4;     %两信号源角度

%第二步： 阵元获取观测数据
A=[1 exp(-j*2*pi*d*sin(xx)/l) exp(-j*2*pi*2*d*sin(xx)/l) exp(-j*2*pi*3*d*sin(xx)/l) exp(-j*2*pi*4*d*sin(xx)/l) exp(-j*2*pi*5*d*sin(xx)/l) exp(-j*2*pi*6*d*sin(xx)/l) exp(-j*2*pi*7*d*sin(xx)/l);
   1 exp(-j*2*pi*d*sin(yy)/l) exp(-j*2*pi*2*d*sin(yy)/l) exp(-j*2*pi*3*d*sin(yy)/l) exp(-j*2*pi*4*d*sin(yy)/l) exp(-j*2*pi*5*d*sin(yy)/l) exp(-j*2*pi*6*d*sin(yy)/l) exp(-j*2*pi*7*d*sin(yy)/l); 
    ]';
S=[sin(2*pi/7*[0:N_s-1]);sin(w*[0:N_s-1])];      %2行N_s列矩阵。!!注意！！起始值从0开始的，所以为了矩阵长度一致，到N_s-1结束
%x=A*S+randn(M,N_s)+j*randn(M,N_s);                %%M行N_s列矩阵，观测数据的获得，加的是白噪声。 第一种高斯白噪声方法，通过randn()，同时需要考虑信噪比snr.;
x1=A*S;
x=awgn(x1,SNR);                                    %                                           第二种加高斯白噪声的方法，通过awgn();

%第三步：对x取自相关，并对其特征分解
Rx=x*x';                                           %根据自相关函数的定义，求自相关 !!注意：x'与conj(x)是两个不同的概念。
[V D]=eig(Rx);
[B index]=sort(diag(D));                           %对特征值进行排序， eig默认情况下是升序排列的。换言之，中的特征值是从小到大升序。本语句是在一维下进行的sort
En=V(:,index(1:M-N));                              %提取噪声空间

%第四步：空间域进行扫描
Fs=100;                                            %要满足抽样定理。因为信号的最高频率为10HZ,由抽样定理可取抽样频率的取值
step=1/Fs;                                         %步进值和抽样频率存在倒数关系。
theta=-90:step:90;                                 %在【-pi/2，pi/2】范围内进行扫描
for k=1:length(theta)                                   
    AA=[1 exp(-j*2*pi*d*sin(theta(k)*pi/180)/l) exp(-j*2*pi*2*d*sin(theta(k)*pi/180)/l) exp(-j*2*pi*3*d*sin(theta(k)*pi/180)/l) exp(-j*2*pi*4*d*sin(theta(k)*pi/180)/l) exp(-j*2*pi*5*d*sin(theta(k)*pi/180)/l) exp(-j*2*pi*6*d*sin(theta(k)*pi/180)/l) exp(-j*2*pi*7*d*sin(theta(k)*pi/180)/l)]';
                                                   %此处A为8行1列矩阵。但根据谱的公式中，可对上述中的AA进行技巧性的设置。                                                                           
    P=AA'*En*En'*AA;
    Pmusic(k)=abs(1/P);                            %此处的PPmusic(k)不能用PPmusic代替，因为这是在循环句中，若被替换后，PPmusic的值会不断覆盖，直到最后一次循环，而此时的PPmusic只是一个值                                                             
end
Pmusic=20*log10(Pmusic);                            %转换为dB单位
sita=-90:step:90;
plot(sita,Pmusic); 
grid on
title('基于二阶距的music算法');
xlabel('角度');
toc;