tic;
clear all
j=sqrt(-1);
N=21;
SNR_Ori=-13:11;
fs=4800;
%xiuzheng=[1 1/1.5819]; % for broadside seven
xiuzheng=[(1.5445*1.0766*1.5849) 1.1523]; % for broadside 21
for aa=1:1
for g=1:length(SNR_Ori)
% [aa g length(SNR_Ori)] % 用于观察程序的进度
SNR(g)=10^(SNR_Ori(g)/10);
f=2000;
c=1500;
k=2*pi*f/c;
lmda=c/f;
theta0=0.4*pi;
t=0:1/fs:2/f;
Ls=length(t);
% D=lmda*[-1.5 -1 -0.5 0 0.5 1 1.5;
% -1.5 -1.34 -1.17 0 1.17 1.34 1.5];
D=lmda*[-5 -4.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5;
-5 -4.82 -4.49 -4.12 -3.74 -3.34 -2.96 -2.61 -2.36 -0.2 0 0.2 2.36 2.61 2.96 3.34 3.74 4.12 4.49 4.82 5];
for p=1:N
for q=1:N
dd(p,q)=besselj(0,k*abs(D(aa,p)-D(aa,q)));% dd 代表 不同阵元之间的间距。
end
end
W=eye(N,N);
Q=(1000/1001)*dd+(1/1001)*W; % Q
v=zeros(N,1);
v_diff=zeros(N,1);
for n=1:N
v(n)=exp(j*k*D(aa,n)*sin(theta0));
v_diff(n)=exp(j*k*D(aa,n)*sin(theta0))*(j*k*D(aa,n)*cos(theta0));
end
v_white=Q^(-0.5)*v;
v_projection=W-v_white*(v_white'*v_white)^(-1)*v_white';
v_white_diff=Q^(-0.5)*v_diff;% 这是对白化后的导引向量求导数
fanshu1=norm(v_projection*v_white_diff);
fanshu2_square=norm(v_white)^(2);
FI=2*SNR(g)*fanshu1^(2)*SNR(g)*fanshu2_square/(1+SNR(g)*fanshu2_square);
CRB(aa,g)=10*log10(1/(FI*Ls));
var_noise=1;
for n=1:Ls
arfa(n,1)=conj(exp(j*2*pi*f*(n-1)/fs));
end
for b=1:1000%蒙特卡罗循环开始
for n=1:N
noise_white(n,:)=wgn(1,Ls,((1/1001)*var_noise),'linear');
end
MU=zeros(1,N);
SIGMA=(1000/1001)*var_noise*dd;
noise_iso= mvnrnd(MU,SIGMA,Ls);
noise_iso=noise_iso.';
noise=noise_white+noise_iso; %每次循环都可以重新产生一次噪声
noiseP=noise*noise'/Ls;
for gg=1:N
var_signal(gg,g)=SNR(g)/xiuzheng(aa)*noiseP(gg,gg);
end
theta_change=-0.4*pi:0.001*pi:0.4*pi;
for m=1:length(theta_change)
for n=1:N
vector_change(n,1)=exp(j*2*pi*D(aa,n)*sin(theta_change(m))/lmda); %这还要改变阵
end
for a=1:Ls
for gg=1:N
Sig(gg,a)=sqrt(var_signal(gg,g))*arfa(a)*v(gg);% 记得这要变 var_signal
end
end
Xs=Sig+noise;
for n=1:Ls
L(n)=(norm(Xs(:,n)-vector_change.*sqrt(var_signal(:,g))*arfa(n)))^(2); %s_amplitude 要随信噪比变
end
LL(m)=sum(L);% 这里是需要修改的地方
end
theta_estimation(b)=theta_change(find(LL==min(LL)));
end
mse_ori(aa,g)=sum((theta_estimation-theta0).^(2))/1000;
mse(aa,g)=10*log10(mse_ori(aa,g));
end
end
% plot(SNR_Ori,CRB(1,:),'--r',SNR_Ori,CRB(2,:),'k',SNR_Ori,mse(1,:),'ro',SNR_Ori,mse(2,:),'kd');
% axis([-13 11 -60 10]);
% xlabel('Input SNR(dB)');
% ylabel('bearing CRB and MSE(dB ref rad^2)');
% legend('---CRB(UA)','─ CRB(INOA)','○ MSE(UA)','◇ MSE(INOA)');
% plot(SNR_Ori,CRB_ua,'--r',SNR_Ori,CRB_inoa,'k',SNR_Ori,mse_ua,'ro',SNR_Ori,mse_inoa,'kd');
% axis([-13 11 -60 10]);
toc;
评论0