function X = CRLB( BSN, MSP, R )
%   Cramer-Rao Lower Bound 是無偏估計的理論下界
%   CRML參數說明:
%       BSN: 基站的個數；
%       MS:  移動臺的位置，其中MSx/MSy均在[0，1]；
%       R:   小區半徑。
%
%       X:   CRLB是一矩陣的形式，主對角線上第一個元素表示橫坐標的CRLB；
%            第二個元素表示定位精度的CRLB。
%   See also CRLB.m

%   李金倫，西南交通大學
%       10 December, 2004, 第一版

%   參數檢查：
if  nargout>1,
    error('Too many output arguments.');
end
if nargin<2 | nargin>3,
    error('Wrong number of input arguments.');
end
% if BSN <= 3,
%     error('The number of BSs must be larger than 3 for this program.');
% end
flag = size(MSP);
if flag(1)~=1 | flag(2)~=2,
    error('Wrong position vector!');
end

%   初始參數設置:
%     BSN = 4;
%     MS  = [0.5, 0.8];
%     R   = 5000;
    BS  = [0, sqrt(3)*R, 0.5*sqrt(3)*R, -0.5*sqrt(3)*R, -sqrt(3)*R, -0.5*sqrt(3)*R, 0.5*sqrt(3)*R;
           0, 0,         1.5*R,         1.5*R,          0,          -1.5*R,         1.5*R];
    c   = 3*10^8;
    Dev = 900/(c*c);
    Q   = 0.5*Dev*(eye(BSN-1) + ones(BSN-1));
    %Q = Dev*eye(BSN-1);
%   算法流程：
    MS = R*MSP;
    
    BSR_1  = sqrt(MS(1)*MS(1) + MS(2)*MS(2)); % R(1)
    
    % R(2) --- R(BSN)
    for i = 1: BSN-1,
        BSR(i) = sqrt((BS(1, i+1) - MS(1))*(BS(1, i+1) - MS(1)) + (BS(2, i+1) - MS(2))*(BS(2, i+1) - MS(2)));
    end
%    % d(2) --- d(BSN)
%     for i = 1: BSN-1,
%         d(i) = (BSR(i) - BSR_1)/c;
%     end
%     Sigma = cov(d)
%     Q = 0.5*Sigma*(eye(BSN-1) + ones(BSN-1));
    % Ga0
    for i = 1:BSN -1,
        Ga(i,1) = -BS(1, i+1);
        Ga(i,2) = -BS(2, i+1);
        Ga(i,3) = -BSR(i);
    end
    % Ga'
    mGa = [1, 0; 0, 1; 1, 1];
    % B
    B = zeros(BSN-1, BSN-1);
    for i = 1: BSN-1,
        B(i, i) = BSR(i);
    end
    % B'
    mB = [MS(1), 0, 0; 0, MS(2), 0; 0, 0, BSR_1];
    % B''
    mmB = [MS(1), 0; 0, MS(2);];
    
%   輸出
    Crlb = c*c*inv(mmB*transpose(mGa)*inv(mB)*transpose(Ga)*inv(B)*inv(Q)*inv(B)*Ga*inv(mB)*mGa*mmB);
    
    if nargout == 1, 
        X = Crlb;
    elseif nargout == 0,
        disp( Crlb );
    end