statistics_classical_music_zhenyuan.m

含有波束空间MUSIC算法、阵元数变化的经典MUSIC算法统计分析、与平滑相关的矩阵分解MUSIC算法

%MUSIC ALOGRITHM
%DOA ESTIMATION BY CLASSICAL_MUSIC
clear all;
close all;
clc;

bbb=zeros(3,10);
for kk=1:10
sensor_number=[3 4 6 8 10 12 14 16 18 20];%阵元数


aaa=zeros(1,300);
for k=1:300


source_number=1;%信元数

N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=pi/4;%信号频率
l=2*pi*3e8/w;%信号波长  
d=0.5*l;%阵元间距
snr=-10;%信噪比(dB)

source_doa=50;%信号的入射角度
A=[exp(-j*(0:sensor_number(kk)-1)*d*2*pi*sin(source_doa*pi/180)/l)].';

s=(10.^(-1/2))*exp(j*w*[0:N_x-1]);%仿真信号
%x=awgn(s,snr);
x=A*s+(1/sqrt(2))*(randn(sensor_number(kk),N_x)+j*randn(sensor_number(kk),N_x));%加了高斯白噪声后的阵列接收信号

R=x*x'/snapshot_number;

%[V,D]=eig(R);
%Un=V(:,1:sensor_number(kk)-source_number);
%Gn=Un*Un';
[V,D]=eig(R);
D=diag(D);
disp(D);
Un=V(:,1:sensor_number(kk)-source_number);
Gn=Un*Un';

searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度
 for i=1:length(searching_doa)
   a_theta=exp(-j*(0:sensor_number(kk)-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l);
   Pmusic(i)=1./abs((a_theta)'*Gn*a_theta);
 end
plot(searching_doa,10*log(Pmusic));
%axis([-90 90 -90 90]);
xlabel('入射角/度');
ylabel('谱峰');
legend('Music Spectrum');
title('经典MUSIC估计');
grid on;
%开始统计分析,先找出谱峰所对应的极值点所对应的估计信号入射角度
aa=diff(Pmusic);
aa=sign(aa);
aa=diff(aa);
bb=find(aa==-2)+1;
[t1 t2]=max(Pmusic(bb));
estimated_source_doa=searching_doa(bb(t2));
disp('estimated_source_doa');
disp(estimated_source_doa);

aaa(:,k)=estimated_source_doa;

disp('estimated_source_doa');
disp(estimated_source_doa);
end
disp(aaa);

%求解均方根误差和标准偏差
E_source_doa=sum(aaa(1,:))/300;%做300次试验的均值
disp('E_source_doa');
disp(E_source_doa);
E_source_doa=mean(aaa);
disp(E_source_doa);
MSE_source_doa=sum((aaa(1,:)-source_doa).^2)/300;%做300次试验的均方误差
disp('MSE_source_doa');
disp(MSE_source_doa);
var_source_doa=sum((aaa(1,:)-E_source_doa).^2)/300;%做300次试验的方差
disp('var_source_doa');
disp(var_source_doa);
deviation_source_doa=E_source_doa-source_doa;%做300次试验的标准偏差
disp('deviation_source_doa');
disp(deviation_source_doa);

bbb(:,kk)=[deviation_source_doa,var_source_doa,MSE_source_doa].';
end
disp(bbb);

figure(2);

plot(sensor_number,bbb(3,:),'r');
hold on
plot(sensor_number,bbb(2,:),'b');
hold on
plot(sensor_number,bbb(1,:),'k');

title('经典MUSIC统计分析');
xlabel('阵元数');
grid on;