haoren.m

谱估计源程序：是谱估计中的经典算法基于线阵MUSIC等算法的源程序

%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');