%classical high resolution DOA estimation algorithm
%including Bartlett beamforming algorithm;Capon algorithm;Esprit;
%MUSIC; WSF; etc;
clear all; close all;clc;
J=sqrt(-1);
source_number=3;
source_doa=[30 40 50];
sensor_number=8;
snapshot_number=100;
snr=10;
A=exp(-J*(0:sensor_number-1)'*pi*sin(source_doa*pi/180));
s=(randn(source_number,snapshot_number)+J*randn(source_number,snapshot_number))/sqrt(2);
x=A*s;
y=awgn(x,snr);
R=y*y'/snapshot_number;
[V,D]=eig(R);
Un=V(:,1:sensor_number-source_number);
Gn=Un*Un';
searching_doa=0:0.1:90;
for i=1:length(searching_doa)
a_theta=exp(-J*(0:sensor_number-1)'*pi*sin(pi*searching_doa(i)/180));
P_BF(i)=abs((a_theta'*R*a_theta)./(a_theta'*a_theta));
P_capon(i)=1./abs((a_theta'*inv(R)*a_theta));
P_music(i)=1./abs((a_theta'*Gn*a_theta));
end
figure(1);
plot(searching_doa,P_BF);
legend('Bartlett spectrum');
figure(2);
plot(searching_doa,P_capon);
legend('Capon spectrum');
figure(3);
plot(searching_doa,P_music);
legend('Music spectrum');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
程序2:以下是分布式目标参数估计算法
%Array signal processing
%The new method to estimate DOA of Distributed signals
%The proposed iteration method
%2005.9.21 Guoxiansheng
clc; clear all;close all;
warning off MATLAB:divideByZero
NN=100;%The number of snapshots
NM=32;% The number of array units
snm=2;%Signal numbers
J=sqrt(-1);%Complex number unit
ag=zeros(snm,1);%The DOA of the signals
ag(1)=-2.5; ag(2)=-1.5;ag(3)=-2.5;
delta1(1)=2;delta1(2)=10;delta1(3)=3;
rho(1)=0.9;rho(2)=0.8;rho(3)=0.7;
%delta(1)=0;delta(2)=0;delta(3)=0;
agg=(ag/180)*pi;
delta=(delta1/180)*pi;
SNR=zeros(1,snm);%SNR of the snm signals
for i=1:snm
SNR(i)=10;
end
Pn=10^-4; %noise power
An=sqrt(Pn/2); %noise amplitude
Ps=10.^(SNR/10)*Pn;% signal power
Amp=sqrt(Ps);%signal amplitude
AA=zeros(NM,snm);%The array manifold
AA(1,:)=1;
noise=zeros(NM,1);
for i=1:snm
for j=2:NM
%AA(j,i)=AA(j,i)+rho(i)^(j-1)*exp(J*pi*(j-1)*sin(agg(i)));
%AA(j,i)=AA(j,i)+2*(1-cos(j-1)*delta(i))/((j-1)*delta(i)).^2*exp(J*pi*(j-1)*sin(agg(i)));
AA(j,i)=AA(j,i)+exp(-(delta(i)*(j-1))^2/2)*exp(J*pi*(j-1)*sin(agg(i)));
%AA(j,i)=AA(j,i)+sinc((j-1)*delta(i))*exp(J*pi*(j-1)*sin(agg(i)));
end
end
fig=figure;plot(abs(fftshift(fft(AA))));zoom xon;
sigs=zeros(snm,NN);%To generate the uncorrelated narrow band signals
for i=1:snm
for j=1:NN
sigs(i,j)=Amp(i)*(randn+J*randn);
end
end
yss=zeros(NM,NN);%The signals with noise
for j=1:NN
noise=An*(randn(NM,1)+J*randn(NM,1));
yss(:,j)=AA*sigs(:,j)+noise;
end
%%%%%%%%%%%%%%%%%%%%%%%
%Traditional MUSIC
Rx=zeros(NM,NM);
for i=1:NN
Rx=Rx+yss(:,i)*yss(:,i)';
end
Rx=Rx/NN;
[evector,evalue]=eig(Rx);
evalue1=diag(evalue);
Gn=evector(:,1:NM-snm);
Us=evector(:,NM-snm+1:NM);
Gnn=Gn*Gn';
NNN=160;qq=16;
for i=1:NNN
agv=(-(i-NNN)/180/qq)*pi;
for j=1:NM
Gs(j)=exp(J*pi*(j-1)*sin(agv));
end
pp(i)=1/abs(Gs*Gnn*Gs');
end
pp=pp/max(pp)
fig=figure;plot(((1:NNN)-NNN)/qq,pp);zoom xon;
title('traditional MUSIC');
%%%%%%%%%%%%%%%%%
%Proposed
for i=1:NNN
agv=(-(i-NNN)/180/qq)*pi;
for j=1:NM
Gs(j)=(exp(J*pi*(j-1)*sin(agv)));
end
pp1=inv(real(diag(Gs)*Gnn*diag(Gs')));
[dde,ddv]=eig(pp1);
%pp1=(diag(Gs)*((Gnn))*diag(Gs'));
pp2=diag(ddv);
for m=1:NM
[vvi,vvv]=max(pp2);
pp3(m,i)=log10(vvi);
vv(vvv)=0;
end
end
for m=1:snm
pp3(m,:)=-(-1)^m*pp3(m,:)/max(pp3(m,:));
end
%est_num1=detect_peak(pp3);
%est_doa1=(est_num1-NNN)/qq;
fig=figure;
plot(((1:NNN)-NNN)/qq,pp3,'-d');zoom xon;
legend('pp2(1)',2);
% fig=figure;
% plot(((1:NNN)-NNN)/qq,pp-pp3,'-d');zoom xon;
title('Proposed MUSIC');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%maxmax
for i=1:NNN
agv=((i-NNN)/180/qq)*pi;
for j=1:NM
Gs(j)=(exp(J*pi*(j-1)*sin(agv)))^2;
end
qq1=Us'*diag(Gs)*(conj(Us));
[s,v]=eig(qq1*conj(qq1));
vv=(diag(v));
for m=1:snm
[vvi,vvv]=max(vv);
qq2(m,i)=log10(1/(1-vvi));
vv(vvv)=0;
end
end
for m=1:snm
qq2(m,:)=-(-1)^m*qq2(m,:);%/max(qq2(m,:));
end
fig=figure;
plot(((1:NNN)-NNN)/qq,qq2);zoom xon;
fig=figure;
plot(((1:NNN)-NNN)/qq,sum(qq2),'-d');zoom xon;
disp('End of programme');