doa_mdalg.m

DOA estimation methods

clc;
%信源
%--------------------------------------------
Pd=2000;
Fd=1;
Fs=4*Fd;
R=0.5;
Delay=5;
M=4;
x1=randint(Pd,1,M); 
x2=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);
y2=modmap(x2,Fd,Fs,'qask',M);
[rcv_a1,ti]=rcosflt(y1,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a2,ti]=rcosflt(y2,Fd,Fs,'fir/sqrt/Fs',R,Delay);

s1=amodce(rcv_a1,10,'qam');
% s2=amodce(rcv_a2,10,'qam');
s2=s1; %假设信号s1和s2相干
save sig3 s1 s2
%--------------------------------------------

%全局变量定义
%--------------------------------------------
clear
m=8; %阵元数
angle1=-10;
angle2=10;
th=[angle1;angle2];
nn=1024;
SN1=10;
SN2=10;
SN3=12;
sn=[SN1;SN2];
degrad=pi/180;
%--------------------------------------------

%生成离散信号
%--------------------------------------------
load sig3 
tt=1:nn;
S=[s1(tt).';s2(tt).'];
nr=randn(m,nn);
ni=randn(m,nn);
U=nr+j*ni; %每路噪声均值为2.
Ps=S*S'/nn;
ps=diag(Ps);
refp=2*10.^(sn/10);
tmp=sqrt(refp./ps);
S2=diag(tmp)*S;
%--------------------------------------------

%协方差矩阵及前后向平滑
%--------------------------------------------
tmp2=[0:m-1]'; %阵元位置
A=exp(-i*pi*tmp2*sin(th'*degrad)); %方向矩阵
X=A*S2+U;
Rxx=X*X'/nn;
a=eye(8);
b=a(:,8:-1:1);%构造交换矩阵
RXX=Rxx-b*Rxx.'*b;
%--------------------------------------------

%奇异值分解
%--------------------------------------------
[Q,SS,W]=svd(RXX); 
Vu=Q(:,3:8); %小特征值对应的特征向量，即噪声子空间
%--------------------------------------------

%做图
%--------------------------------------------
th2=[-90:90]';
tmp2=[0:m-1]';
a2=-i*tmp2*pi*sin(th2'*degrad);
A2=exp(a2);
num=diag(A2'*A2);
Ena=Vu'*A2;
den=diag(Ena'*Ena);
doa=num./den;
semilogy(th2,abs(doa));
axis([-90 90 0.1 1e5]);
grid on;
hold on;
%--------------------------------------------
clear all;
pi=3.1415926;% pi
fs=2.5e+12;% modulation rate
tp=1e-5;% pulse duration
f0=0.1e+8;% shift frequence 
f1=0.3e+8;
f2=0.5e+8;
f3=0.7e+8;
BW=2*tp*fs+f0-f0;% bandwidth
tal=1/2/BW;% sampling interval
t=0:tal:tp;% time duration of pulse
%t1=tp*2:tal:tp*3;
ns=round(tp/tal)+1;% sampling number of transmitted signals
deltf=BW/50;% frequence interval for FT computaion
nf=60;% FT point number
sts0=exp(j*2*pi*f0*t+j*2*pi*t.^2*fs);
sts1=1*exp(j*2*pi*f1*t+j*2*pi*t.^2*fs);
sts2=1*exp(j*2*pi*f2*t+j*2*pi*t.^2*fs);
sts3=1*exp(j*2*pi*f3*t+j*2*pi*t.^2*fs);
%sts=sts1+sts2+sts3+sts0;
sts=sts0;% signal model
stsft=zeros(1,nf);% FT of signal
%% FT computation of signal
for m=1:nf
    for n=1:ns
        stsft(m)=stsft(m)+(sts(1,n))*exp(-j*2*pi*(m-0)*deltf*(n)*tal);
    end
    psdsts(m)=stsft(m)*stsft(m)';% PSD computation of signal
end
absstsft=abs(stsft);
fd=(0:nf-1)*deltf;
figure(1)
plot(fd,absstsft);
figure(2)
plot(psdsts);
phast=zeros(1,ns);% phase of transmitted signals
phast=phase(sts);%+rand(1,Numst)*2*pi;
ampst=zeros(1,ns);%amplitude of transmitted signals
ampst=abs(sts);
Rmatrix=zeros(ns,ns);% R2 matrix
Numstopband=1;% stopband number
freqstopband=zeros(Numstopband,2);% stopband frequence 
%freqstopband=[0,9;16,29;36,49;56,69;76,130]*(1e+6);
freqstopband=[21,32]*(1e+6);
%weights
wit=zeros(1,Numstopband);
wit=[0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1];
tsample=tal;
 % R2 matrix
for k=1:ns
    for m=1:ns
        if k==m
            for n=1:Numstopband
                Rmatrix(k,m)=Rmatrix(k,m)+0.1*freqstopband(n,2)-0.1*freqstopband(n,1);
            end
        else
            for n=1:Numstopband
                Rmatrix(k,m)=Rmatrix(k,m)+0.1*(exp(-2*pi*j*(k-m)*freqstopband(n,2)*tsample)-exp(-2*pi*j*(k-m)*freqstopband(n,1)*tsample))/(-2*pi*j*(k-m)*tsample);
            end
        end
    end
end
%diagonal loading and normalizing Rmatrix;
Mload=eye(ns)*Rmatrix(1,1)*1e-5;
Rmatrix=Rmatrix+Mload;
Rmatrix=Rmatrix/Rmatrix(1,1);
penfun=sts*Rmatrix*sts'
% % To find a proper mu
uu=eig(Rmatrix);
ul=1/max(uu);