music.m

波束形成中DOA估计MUSIC算法 matlab程序

M=100; %接收阵元个数
N=100; %快拍次数
end_f=0.5;%搜索截止频率
count=0;
for delta_f =0:0.001:end_f
count=count+1;
end
Pw=zeros(count,20);
%产生数据
for loop=1:20
w = randn(2000,1);
for n=1:300
x(n)=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+w(n+50*loop);
end
%构造X矩阵
X=zeros(M,N);
for Nn=1:N
X(:,Nn)=x(Nn:Nn+M-1)';
end
%求相关矩阵
Rxx=zeros(M,M);
for Nn=1:N
Rxx=Rxx+X(:,Nn)*X(:,Nn)';
end
Rxx=Rxx/N;
%求特征值、特征向量及阶数P
[U,SS,V]=svd(Rxx);
bottom = 0;
for hh=1:M
bottom = bottom + SS(hh,hh)^2;
end
overhead=0;
v = overhead/bottom;
p=0;
while v<0.999
p = p + 1;
overhead = overhead + SS(p,p)^2;
v = overhead/bottom;
end
%
end
p
S=U(:,1:p);
G=U(:,p+1:M);
%构造M阶矩阵 
i=0;
for delta_f =0:0.001:end_f
i=i+1;
f(i)=delta_f;
w(i)=2*pi*delta_f;
for Mm=1:M
aw(Mm)=exp(-j*(Mm-1)*w(i));
end
%用噪声子空间构造功率谱函数
Pw(i,loop)=1/(aw*(G*G')*aw');
end
%求空间谱模值、平均值、绘功率谱图
absPw=20*(log(abs(Pw)));
Pw_mean=mean(abs(Pw),2);
plot(f,Pw_mean);
%找出功率谱的峰值
k=0;
for i=2:count-1
if Pw_mean(i) > Pw_mean(i-1) && Pw_mean(i) > Pw_mean(i+1)
k=k+1;
f_max(k)=f(i);
Pw_max(k)=Pw_mean(i);
end
end
%找出峰值点对应的频率
Pw_max=sort(Pw_max);
f1=0;
f2=0;
for i=1:count
if Pw_mean(i)==Pw_max(k)
f2=f(i);
Pw_mean(i)
end
end
for i=1:count
if Pw_mean(i)==Pw_max(k-1)
f1=f(i);
Pw_mean(i)
end
end
f1
f2