function CRLB = TDOAFDOALocMultiDisjMvgSrcSenPosErrCRLB(so,s_doto,uo,u_doto,Q_beta,Q_alpha)
%
% This program calculates the CRLB of multiple disjoint source positions
% and velocities when the sensor positions and velocities are subject to 
% random noises. The positioning parameters are TDOA and FDOA. 
%
% The derivation of the CRLB can be found in Ming Sun and K. C. Ho, 
% "An Asymptotically Efficient Estimator for TDOA and FDOA Positioning of 
% Multiple Disjoint Sources in the Presence of Sensor Location Uncertainties,"
% IEEE Trans. Signal Process., vol. 59, pp. 3434-3440, Jul. 2011.
%
% Usage: CRLB = TDOAFDOALocMultiDisjMvgSrcSenPosErrCRLB(so,s_doto,uo,u_doto,Q_beta,Q_alpha)
%
% Input Parameter List
% so:  2xM or 3xM true sensor position matrix.
%      M: number of sensors.
%      so(:,i) is the true position of the ith sensor.
% s_doto: 2xM or 3xM true sensor velocity matrix.
%         s_doto(:,i) is the true velocity of the ith sensor.
% uo: 2xK or 3xK true source position. 
%     K: the number of the sources to be localized.
% u_doto: 2xK or 3xK true source velocity.
% Q_beta: 4Mx4M, or 6Mx6M covariance matrix of the sensor position and velocity
%                       reshape([s,s_dot],4M,1) or reshape([s,s_dot],6M,1).
% Q_alpha: 2(M-1)x2(M-1) covariance matrix of the TDOAs and FDOAs.
%
% The program returns the 4Kx4K or 6Kx6K CRLB of [uo;u_doto].
%
% The program can be used for 2D(Dim=2) or 3D(Dim=3) localization.
%
%
% Ming Sun and K. C. Ho      03-05-2011
%
%       Copyright (C) 2011
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

M = size(so,2);           % Number of sensors.
N = size(so,1);           % Dimension of the localization problem.

K = length(uo)/N;         % Number of sources.

E = repmat(eye(N),M,1);

dalpha_dtheta = [];
dalpha_dbeta = [];
for k = 1 : K
    a1 = (uo(N*(k-1)+1:N*k)-so(:,1))'/norm((uo(N*(k-1)+1:N*k)-so(:,1)));
    
    b1 = (u_doto(N*(k-1)+1:N*k)-s_doto(:,1))'/norm(uo(N*(k-1)+1:N*k)-so(:,1)); 
    temp = (u_doto(N*(k-1)+1:N*k)-s_doto(:,1))'*(uo(N*(k-1)+1:N*k)-so(:,1))*(uo(N*(k-1)+1:N*k)-so(:,1))';
    temp = temp/norm(uo(N*(k-1)+1:N*k)-so(:,1))^3;
    b1 = b1 - temp;
    
    for i = 2 : M         % Generate the partial derivatives.
        C1(i-1,:) = [zeros(1,(i-1)*N),(uo(N*(k-1)+1:N*k)-so(:,i))'/norm((uo(N*(k-1)+1:N*k)-so(:,i))),...
            zeros(1,N*(M-i))];
        C1(i-1,1:N) = -a1;
        
        bi = (u_doto(N*(k-1)+1:N*k)-s_doto(:,i))'/norm(uo(N*(k-1)+1:N*k)-so(:,i)); 
        temp = (u_doto(N*(k-1)+1:N*k)-s_doto(:,i))'*(uo(N*(k-1)+1:N*k)-so(:,i))*(uo(N*(k-1)+1:N*k)-so(:,i))';
        temp = temp/norm(uo(N*(k-1)+1:N*k)-so(:,i))^3;
        bi = bi - temp;
        C2(i-1,:) = [zeros(1,(i-1)*N),bi,zeros(1,N*(M-i))];
        C2(i-1,1:N) = -b1;
    end;
    dr_du = [[C1;C2]*E,[zeros(M-1,N);C1*E]];
    dr_ds = [-[C1;C2],[zeros(M-1,M*N);-C1]];
    
    dalpha_dtheta = blkdiag(dalpha_dtheta,dr_du);
    dalpha_dbeta = [dalpha_dbeta;dr_ds];
end;

X = dalpha_dtheta'*inv(Q_alpha)*dalpha_dtheta;
Y = dalpha_dtheta'*inv(Q_alpha)*dalpha_dbeta;
Z = inv(Q_beta) + dalpha_dbeta'*inv(Q_alpha)*dalpha_dbeta;

CRLB = inv([X,Y;Y',Z]); 
CRLB = CRLB(1:2*K*N,1:2*K*N); % Return the CRLB.