www.pudn.com > SSR_DOA.rar > SSR_doa.m, change:2014-03-13,size:3436b

clc;
clear;
close;

lambda=1;
d=lambda/2;      %阵元间距离，取为入射波长的一半
K=500;            %采样快拍数
theta=[-50,40];  %入射角度
SignalNum=length(theta);     %入射信号数量
Nnum=5;                 %%阵列阵元数量
SNR1=-8;          %%信噪比
Aratio=sqrt(10^(SNR1/10));    %信号幅度与噪声幅度比值，并假设信号幅度为1

Fs=5*10^3;               %信号频率
Fc=[2*10^3,5*10^3,8*10^3];       %入射信号频率
fs=20*10^3;

thetatest=(-90*pi/180:1*pi/180:90*pi/180);   %theta角度搜索范围
thetanum=length(thetatest);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算信号协方差矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_Vector=(1:K)/fs;
A=zeros(Nnum,SignalNum);
SignalVector=zeros(SignalNum,K);
%NoiseVector=zeros(Nnum,K);
Xt=zeros(Nnum,K);
%%构造A矩阵
for k2=1:SignalNum
for k1=1:Nnum                          %1:12
At(k1)=exp(j*(k1-1)*2*pi*d*sin(theta(k2)*pi/180)/lambda);
A(k1,k2)=At(k1);
end
end
%%%构造信号矩阵和噪声矩阵
for k1=1:SignalNum
SignalVector(k1,:)=exp(j*2*pi*Fc(k1).*T_Vector);  %信号
end
Xtt=A*SignalVector;

%NoiseVector=sqrt(0.5)*(randn(Nnum,K)+j*randn(Nnum,K));
for kk=1:Nnum
Xt(kk,:)=awgn(Xtt(kk,:),SNR1,'measured');
end
Rx=(Xt*Xt')./K;

Rs=(SignalVector*SignalVector')./K;
sigm_s=Rs(:,1);
% %%%%%%%%%%%%%%%%%%%%%%%%%%-----特征值分解----%M%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[V, D] = eig(Rx);                %   X*V = V*D
DD     = diag(D);               % 对角阵变矢量
% [DD idx] = sort(DD, 'descend'); % 按从大往小排序特征值
Un     = V(:, 1:Nnum-SignalNum);         % 噪声子空间
Us     = V(:, Nnum-SignalNum+1 : end);  % 信号子空间
e1=[1,zeros(1,Nnum-1)].';
sigm_n=min(DD);        %最小特征值^作为 的估计
I=eye(Nnum);

for k1=1:thetanum
Atemp0=exp(j*2*pi*d/lambda*sin(thetatest(k1))*[0:Nnum-1]).';
S(k1)=1/(Atemp0'*Un*Un'*Atemp0);
end
% figure(1)
% plot(thetatest.*180./pi,10*log10(abs(S)/(max(abs(S)))));%输出功率（dB）
% grid on;
% % grid on;
% title('Music')
% xlabel('方位角（度）')
% ylabel('输出功率（dB）')

%%%%%%%%%%%%%%%%%%%%%%%%%%----构造--G--selection矩阵%%%%%%%%%%%%%%

M=Nnum;
G=zeros(M*M,2*M-1);
J0=eye(M);
G(:,M)=J0(:);
%for k=1:M-1
for i=1: M-1
J=[zeros(M-i,i),eye(M-i);zeros(i,i),zeros(i,M-i)];
G(:,M-i)=J(:);
J1=J';
G(:,M+i)=J1(:);
end

%%%%%%----Bthita------%%%%
Bthita=zeros(2*M-1,thetanum);
Bt=zeros(1,2*M-1);
for k2=1:thetanum                    %相当于文章thita1---thitaQ
for k1=1:M
Bt(1,k1+M-1)=exp(-j*(k1-1)*2*pi*d*sin(thetatest(k2))/lambda);
Bt(1,k1)=exp(j*(M-k1)*2*pi*d*sin(thetatest(k2))/lambda);
Bthita(:,k2)=Bt';
end
end
%%%----u---K稀疏矢量----
u=zeros(1,thetanum);
for z=1:SignalNum
u(1,theta(z)+ 90+1)=sigm_s(z);    %应该是等于sigm^2,每个信号的噪声方差？？？？

end
u=u';
%%%%%%%%%%%%%%%%%%%%%%%%%%----cvx运算-------%%%%%%%%%%%%%%%%%%%%%
y=Rx(:);

W12=sqrt(K)*(kron((Rx^(-0.5)).',Rx^(-0.5)));
Q=W12*G*Bthita;
S1=W12*(y-sigm_n*I(:))-Q*u;

beita=sqrt(chi2inv(0.999999,M*M));               %卡方分布M*M
cvx_begin
variable u2(181,1)
minimize( norm(u2,1))
subject to
S=W12*(y-sigm_n*I(:))-Q*u2;
%  S=y-G*Bthita*u2-sigm_n*I(:);
norm(S) <=beita ;
cvx_end
figure()
plot(-90:1:90,u2);

%-----------------------------------------------------------------------