pmusic.m

几种music算法

%-----------------------
%MUSIC ALGORITHM
%DOA Estimation by ULA
clear all; 
close all;
clc
N_x=1024; 
n=20; 
l=1.8;
d=0.5*l;
M=3; 
w=[pi/6 pi/6 pi/6]';
xx=pi/6;yy=pi/6;zz=pi/6; %信号的入射角
B=[1 exp(-j*2*pi*d*sin(xx)/l) exp(-j*2*2*pi*d*sin(xx)/l) exp(-j*2*3*pi*d*sin(xx)/l) exp(-j*2*4*pi*d*sin(xx)/l) exp(-j*2*5*pi*d*sin(xx)/l) exp(-j*2*6*pi*d*sin(xx)/l) exp(-j*2*7*pi*d*sin(xx)/l);
   1 exp(-j*2*pi*d*sin(yy)/l) exp(-j*2*2*pi*d*sin(yy)/l) exp(-j*2*3*pi*d*sin(yy)/l) exp(-j*2*4*pi*d*sin(yy)/l) exp(-j*2*5*pi*d*sin(yy)/l) exp(-j*2*6*pi*d*sin(yy)/l) exp(-j*2*7*pi*d*sin(yy)/l);
   1 exp(-j*2*pi*d*sin(zz)/l) exp(-j*2*2*pi*d*sin(zz)/l) exp(-j*2*3*pi*d*sin(zz)/l) exp(-j*2*4*pi*d*sin(zz)/l) exp(-j*2*5*pi*d*sin(zz)/l) exp(-j*2*6*pi*d*sin(zz)/l) exp(-j*2*7*pi*d*sin(zz)/l)]'; 
xxx=2*exp(j*w*[0:N_x-1]);
x=B*xxx+randn(8,N_x)+j*randn(8,N_x);

R=x*x'/N_x;
[V D]=eig(R);
[lambda,index] = sort((diag(D)));
UU=V(:,index(1:5));
theta=-90:0.1:90;
for i = 1:length(theta)
    AA=[1 exp(-j*2*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*2*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*3*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*4*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*5*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*6*pi*d*sin(theta(i)/180*pi)/l) exp(-j*2*7*pi*d*sin(theta(i)/180*pi)/l)];%方向矢量
    WW=AA*UU*UU'*AA';
    Pmusic(i)=abs(1/WW);
end
Pmusic = 20*log10(Pmusic);
fai=-90:0.1:90;
plot(fai,Pmusic);
grid
xlabel('rad')
%---------------
%部分参数设定
h=60;
w=150;
f=1000;
c=1500;
tao=0.005;
fai0=pi/6;
fai=fai*pi/180;
%传播的单向距离
r=h:1:200*h;
r=(h./(sin(fai)));
%声源级SL(轴向)
SL=170.8+10*log10(n)+10*log10(w);
%散射系数Sb
Sb=-5+10*log10(sin(fai).^2);
%等效束宽tai
y=(asin(0.55*c/(f*n*d)+sin(fai0))-asin(-0.55*c/(f*n*d)+sin(fai0)))*180/(2*pi);
tai=10*log10(y)-12.8;
%---------------
%混响级RL
RL=SL-40*log10(r)+Sb+10*log10(r*c*tao/2)+tai+10*log10((Pmusic.^2)./cos(fai));
%---------------
%画混响级随时间变化曲线RL(t)
figure(2);
plot(2*r/c,RL,'b');
title('波束反向散射信号EL'); 
grid on;