function CRLB = TDOALocSenPosErrCalEmtrCRLB(so,uo,c,Q_alpha,Q_c,Q_beta)
%
% This program calculates the CRLB of a stationary source when the sensor
% positions are subject to random noise and a calibration emitter at
% perfect location is present. The positioning parameter is TDOA.
% The derivation of the CRLB can be found in K. C. Ho and
% Le Yang, "On the Use of a Calibration Emitter for Source Localization in 
% the Presence of Sensor Position Uncertainty," IEEE Trans. Signal
% Processing, Vol. 56, pp. 5758-5772, Dec. 2008.
%
% Usage: CRLB = TDOALocSenPosErrCalEmtrCRLB(so,uo,c,Q_alpha,Q_c,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: 2x1 or 3x1 true source position vector.
% c : 2x1 or 3x1 calibration position vector.
% Q_alpha: (M-1)x(M-1) covariance matrix of the RDs from uo.
% Q_c: (M-1)x(M-1) covariance matrix of the RDs from c.
% Q_beta: 2Mx2M or 3Mx3M covariance matrix of the sensor positions  
%                                       reshape(s,2M,1) or reshape(s,3M,1).
%
% The program returns the 2x2 or 3x3 CRLB of the source location.
%
%
% The program can be used for 2D(Dim=2) or 3D(Dim=3) localization.
%
% Le Yang, K. C. Ho      08-01-2009
%                        10-01-2010, revised
%
%       Copyright (C) 2009
%       Computational Intelligence Signal Processing Laboratory
%       University of Missouri
%       Columbia, MO 65211, USA.
%       hod@missouri.edu
%

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

for j = 2 : M         % Generate the partial derivatives.
    dr_du(j-1,:) = (uo-so(:,j))'/norm(uo-so(:,j))-(uo-so(:,1))'/norm(uo-so(:,1));
    dr_ds(j-1,:) = [(uo-so(:,1))'/norm(uo-so(:,1)),zeros(1,(j-2)*N),-(uo-so(:,j))'/norm(uo-so(:,j)),zeros(1,(M-j)*N)];
    drt_ds(j-1,:) = [(c-so(:,1))'/norm(c-so(:,1)), zeros(1,(j-2)*N),-(c-so(:,j))'/norm(c-so(:,j)),  zeros(1,(M-j)*N)];
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 + drt_ds'*inv(Q_c)*drt_ds;

CRLB = inv([X,Y;Y',Z]); 
CRLB = CRLB(1:N,1:N); % Return the CRLB of the source position.