function CRLB = TDOALocSenPosErrMultiSrcDisjCRLB(so,uo,Q_alpha,Q_beta)
% This program calculates the CRLB of multiple disjoint sources when the 
% sensor positions are subject to random noise. The positioning parameter is TDOA.
%
% The derivation of the CRLB can be found in Le Yang and K. C. Ho, 
% "An approximately efficient TDOA localization algorithm in closed-form 
% for locating multiple disjoint sources with erroneous sensor positions," 
% IEEE Trans. Signal Processing, Vol. 57, pp. 4598-4615, Dec. 2009.
%
% Usage: CRLB = TDOALocSenPosErrMultiSrcDisjCRLB(so,uo,Q_alpha,Q_beta)
%
% Input parameter list:
% so: 2xM or 3xM true sensor position matrix. 
%     M: the number of sensors.
%     so(:,i) is the true position of the ith sensor.
% uo: 2Nx1 or 3Nx1 true source position vector.
%     N: the number of sources to be localized.
% Q_alpha: N(M-1)xN(M-1) covariance matrix of the TDOAs (range differences) from uo.
% Q_beta: 2Mx2M or 3Mx3M covariance matrix of the sensor positions  
%                                       reshape(s,2M,1) or reshape(s,3M,1).
%
% The program returns the 2Nx2N or 3Nx3N CRLB of the source position vector.
%
% The program can be used for 2D(Dim=2) or 3D(Dim=3) localization.
%
%
% Le Yang, K. C. Ho      10-01-2010
%
%       Copyright (C) 2010
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

M = size(so,2);       % Number of sensors.           
D = size(so,1);       % Dimension of the localization problem.
N = length(uo)/D;     % Number of sources to be localized.

for i = 1 : N
    u = uo((i-1)*D+1:i*D);
    for j = 2 : M    
        dr_du((i-1)*(M-1)+(j-1),:) = [zeros(1,(i-1)*D),(u-so(:,j))'/norm(u-so(:,j))-(u-so(:,1))'/norm(u-so(:,1)), zeros(1,(N-i)*D)];
        dr_ds((i-1)*(M-1)+(j-1),:) = [(u-so(:,1))'/norm(u-so(:,1)),zeros(1,(j-2)*D),-(u-so(:,j))'/norm(u-so(:,j)),zeros(1,(M-j)*D)];
    end;
end;

X = dr_du' * inv(Q_alpha) * dr_du;
Y = dr_du' * inv(Q_alpha) * dr_ds;
Z = inv(Q_beta) + dr_ds'* inv(Q_alpha) * dr_ds;

CRLB = inv(X) + inv(X)*Y*inv(Z-Y'*inv(X)*Y)*Y'*inv(X);