%opt2d.m: 全自由度空時自適應處理
%--------------------------------------------------------------------------
%start         : 2004.11.04  AM 10:21   yunrisheng
%Latest change : 
%--------------------------------------------------------------------------

clear
tic
%雜波仿真參數
N = 12;                        % 陣元個數
M = 10;                        % 相干脈沖數
CNR = 30;                      % 雜噪比
beta = 1;                      % 雜波折疊系數(beta = 2*v*T/d)
sita_a = -90:.9:90;         % 雜波單元個數               
sita = sita_a*pi/180;
[NN N_bin] = size(sita);
%目標參數
sita_t = -25;                  % 目標DOA
omiga_t = 0.4;                 % 目標Doppler
SNR = 0;                       % 信噪比

%空間導向矢量和時間導向矢量
%空間頻率和Dopple頻率滿足 omiga_d = beta * omiga_s
omiga_s = pi*sin(sita);      
omiga_d = beta*omiga_s;       

aN = zeros(N,N_bin);
bN = zeros(M,N_bin);

aN = exp(-j*[0:N-1]'*omiga_s)./sqrt(N);
bN = exp(-j*[0:M-1]'*omiga_d)./sqrt(M);


%目標空時信號
aN_t = zeros(N,1);
bN_t = zeros(M,1);

aN_t = exp(-j*pi*[0:N-1]'*sin(sita_t*pi/180))/sqrt(N);
bN_t = exp(-j*pi*[0:M-1]'*omiga_t)/sqrt(M);

S_t = zeros(M*N,1);
S_t = kron(aN_t,bN_t);

%計算雜波協方差矩陣
R = zeros(M*N,M*N);                     
S = zeros(M*N,N_bin);                   
ksai = 10^(CNR/10)*(randn(1,N_bin)+j*randn(1,N_bin))/sqrt(2);               %服從正態分布的隨機幅值，方差為1
for ii = 1:N_bin
    S(:,ii) = kron(aN(:,ii),bN(:,ii));  
    R = R + ksai(ii).*(S(:,ii)*S(:,ii)');       
end

%干擾協方差矩陣，雜噪比為30dB
R = R +eye(M*N);     %CNR = 30dB
inv_R = inv(R);                   %逆矩陣
%求特征值譜
[u s v] = svd(R);                       %特征值分解
figure(1);
plot(10*log10(diag(s)));
title('陣元數N=12, 相干脈沖數M=10');
axis([0 120 -10 50]);
xlabel('特征值數目');
ylabel('特征值(dB)');
grid on

P_f = zeros(N_bin,N_bin);
P_min_var = zeros(N_bin,N_bin);
%求雜波譜
for ii = 1:N_bin
    for jj = 1:N_bin
            SS = kron(aN(:,ii),bN(:,jj));
            P_f(ii,jj) = SS'*R*SS;        %傅氏譜
            P_min_var(ii,jj) = 1./(SS'*inv_R*SS);
    end
end        

%最小方差功率譜
figure(2)
mesh(sin(sita),omiga_d/pi,20*log10(abs(P_min_var)));
title('陣元數N=12, 相干脈沖數M=10');
xlabel('方位余弦');
ylabel('歸一化Dopple頻率');
zlabel('功率(dB)');
grid on
%空時最優權向量
w_opt = inv(R)*S_t./(S_t'*inv_R*S_t);   
%w_opt = inv(RR)*a_t;

%求最優空時響應
for ii = 1:N_bin
    for jj = 1:N_bin
        SSS = kron(aN(:,ii),bN(:,jj));
        res_opt(ii,jj) = SSS'*w_opt;
    end
end
                      
figure(3)
%[X,Y]=meshgrid(omiga_d/pi,sita_a);
mesh(omiga_d/pi,omiga_d/pi,10*log10(abs(res_opt).^2))
title('陣元數N=12, 相干脈沖數M=10');
xlabel('歸一化Dopple頻率');
ylabel('方位余弦');
zlabel('功率(dB)');
grid on
%求最優改善因子
for ii = 1:N_bin
    for jj = 1:N_bin
        SS = kron(aN(:,ii),bN(:,jj));
 %       IF(ii,jj) = SS'*inv_R*SS.*trace(R)./(SS'*SS);
        IF(ii,jj) = SS'*inv_R*SS./(SS'*SS);
    end
end
figure(4)
axis([-1 1 10 55]);
plot(omiga_d/pi,10*log10(abs(IF(101,:))));
xlabel('歸一化Dopple頻率');
ylabel('改善因子(dB)');
grid on

toc