%%%兩個信號均存在時，空域最大熵方法求得信號到達方向的估計值
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)])