cccvw.m

采用前后向平滑以及最小冗余线阵的宽带源DOA估计仿真

%加上冗余线阵
clear all
clc
N=1024;
d=2;
M=4;
a=[70,100];
f=1/16;
h=[0,1,3,4];

%SOI signal
	for k=1:d
      s1=randsn1(N);   
   for n=1:N
      ss(k,n)=s1(n)*cos(2*pi*f*n);          %BPSK    
   end
end

SS=zeros(M,N);
for k=1:d
   S{k}=zeros(M,N);
     mov=zeros(M,1);
      for m=1:M
          mov(m,1)=floor(h(m)/(2*f)*cos(a(k)*pi/180));
      end
    
   for i=1:M
     if mov(i,1)>=0 
        for n=1:N-mov(i,1) 
            S{k}(i,n)=ss(k,n+mov(i,1));
        end
     else 
        for n=-mov(i,1)+1:N
             S{k}(i,n)=ss(k,n+mov(i,1));
        end
     end 
   end
SS=SS+S{k};   
end
W=normrnd(0,0.1,M,N);
y=SS+W;

b=2*f;
R=zeros(9,1);
t=1;
     for p=1:M
         for q=1:M
             i=h(p)+h(q)+1;
           for n=1:N-t
            R(i,t)=R(i,t)+y(p,n)*(y(q,(n+t)))*exp(-j*2*pi*b*(n+t/2));
           end
        R(i,1)=R(i,1)/N;
         end
     end 
 
 Y=zeros(5,5);
 for i=1:5
     for j1=1:5
      Y(i,j1)=R(i+j1-1,1); 
     end
 end
 
 
 [U,V,uu]=svd(Y);
 %取噪声子空间特征向量
 abs(V)
Un=U(:,[d+1:M]);
%谱峰搜索
    dtheta=0.1;   %搜索步长
     n=[1:dtheta:180];
for i=1:length(n),
     aa=zeros(M,1);
      for m=1:5
      aa(m,1)=exp(j*(m-1)*pi*cos(n(i)*pi/180));
      end      
    P(i)=-10*log(abs(aa'*Un*Un'*aa));
   %P(i)=1/(abs(aa'*Un*Un'*aa));
end
i=1:length(n);
figure(1)
plot(n(i),P(i)),grid on;
xlabel('来波方向角')
ylabel('P(dB)')

zoom on
%end