最大熵法doa估计.m

DOA估计（雷达来波方向估计）是雷达的重要理论

%%%两个信号均存在时，空域最大熵方法求得信号到达方向的估计值
L  = 8;                                             %%%%%%%%%%% 8阵元
dLamtaR = 0.5;
SNR = input ('SNR = --->');             %%%%%%%%% 输入信噪比
M = input ('M = --->');                    %%%%%%%%% 用作预测的阵元数目
FI =[-30 30]/180*pi;
THETA = 2*pi*dLamtaR*sin(FI);

k = (0:L-1)';
y = exp(j*k*THETA(1))+exp(j*k*THETA(2));
noise = randn(L,1)+j*randn(L,1);
noise = noise/std(noise);
y = y/std(y)+noise/10^(SNR/20);
y_asterisk = conj(y);

Y1 = zeros(L-M,M);
Y2 = zeros(L-M,M);
for k = 1:L-M
    Y1(k,:) = y(M-1+k:-1:k).';
    Y2(k,:) = y_asterisk(L-M-k+2:L-k+1).';
end
Y = [Y1; Y2];
y0 = [y(M+1 : L); y_asterisk(L-M:-1:1)];

h_FBLP = inv(Y' * Y) * Y' * y0;

g_FBLP = [1; -h_FBLP];
fi = -pi/2 : pi/1440 :pi/2;
zi = exp(j * 2 * pi * dLamtaR * sin(fi));

A_FBLP = polyval(g_FBLP, zi);
P_FBLP = 1 ./(abs(A_FBLP) .* abs(A_FBLP)); 
 
plot(fi*(180/pi), P_FBLP);
title(['FBLP, M = ' int2str(M) '  SNR = ' num2str(SNR)])

figure(2)
polar(fi,P_FBLP);
title(['FBLP, M = ' int2str(M) '  SNR = ' num2str(SNR)])