function D_TDOA_new = Refine_TDOA_distance_matrix(D_TDOA_old, dimension, iteration)
    % Input:
    %   D_TDOA_old: M x (N-1) matrix, where M is the number of sources, 
    %               N is the number of sensors, and its elements are 
    %               \tau_{m,n} given by equation (1) in the paper.
    %   dimension:  the dimension of the positions of the sources and sensors.
    %               dimension = 2 is 2_D and dimension = 3 is 3-D.
    %   iteration:  the number of repititions in step 3 of Algorithm 1
    %               in the paper. The default value of "iteration" is 20.
    % Output:
    %   D_TDOA_new: The refinement of D_TDOA_old.
    %
    % Trung-Kien Le, Dominic K. C. Ho   07-15-2018
    %
    % Trung-Kien Le, K. C. Ho and Trung-Hieu Le, "Rank Properties for Matrices 
    % Constructed from Time Differences of Arrival, " IEEE Trans. Signal
    % Process., vol. 66, No. 13, pp. 3491-3503, Jul. 2018.
    %
    %       Copyright (C) 2018
    %       "Rank Properties for Matrices Constructed from Time Differences of Arrival"
    %       University of Missouri
    %       Columbia, MO 65211, USA.
    %       lekt@missouri.edu  and  hod@missouri.edu
    %
    
    [M,N]        = size(D_TDOA_old); N = N+1;
    
    %===== Step 1 of Algorithm ============================================
    Gamma_global = [D_TDOA_old ; D_TDOA_old(2:M,:).^2 - repmat(D_TDOA_old(1,:).^2,M-1,1)];
    
    %===== Step 2 of Algorithm: Refinement of TDOAs based on Theorem 2 ====
    [U,S,V]      = svd(Gamma_global);
    rank         = min([M+dimension , N-1]);
    Gamma_global = U(:,1:rank)*S(1:rank,1:rank)*V(:,1:rank)';
    D_TDOA_new   = Gamma_global(1:M,:);
    D_TDOA       = D_TDOA_new;
    
    %===== Step 3 of Algorithm: Refinement of TDOAs based on Theorem 1 ====
    if dimension == 2   % 2-D case     
        rank1 = min([6 , N-1]);
        for k = 1:iteration
            cnt=zeros(M,1);
            KK = zeros(M,1000,N-1);
             for m1 = 1:M-3
                for m2 = m1+1:M-2
                    for m3 = m2+1:M-1
                        for m4 = m3+1:M
                            Gamma_local = [D_TDOA([m1,m2,m3,m4],:) ; D_TDOA([m2,m3,m4],:).^2 - repmat(D_TDOA(m1,:).^2,3,1)];
                            [U,S,V]     = svd(Gamma_local); 
                            Gamma_local = U(:,1:rank1)*S(1:rank1,1:rank1)*V(:,1:rank1)';
                            cnt([m1,m2,m3,m4])= cnt([m1,m2,m3,m4])+1;
                            KK(m1,cnt(m1),:) = Gamma_local(1,:);
                            KK(m2,cnt(m2),:) = Gamma_local(2,:);
                            KK(m3,cnt(m3),:) = Gamma_local(3,:);
                            KK(m4,cnt(m4),:) = Gamma_local(4,:);
                        end
                    end
                end
            end
            KK = KK(:,1:cnt(1),:);
            D_TDOA = squeeze(mean(KK,2));
        end;
        D_TDOA_new = D_TDOA;
    elseif dimension == 3  % 3-D case       
         rank1 = min([8 , N-1]);
         for k = 1:iteration   
            cnt=zeros(M,1);
            KK = zeros(M,1000,N-1);
            for m1 = 1:M-4
                for m2 = m1+1:M-3
                    for m3 = m2+1:M-2
                        for m4 = m3+1:M-1
                            for m5 = m4+1:M
                                Gamma_local = [D_TDOA([m1,m2,m3,m4,m5],:) ; D_TDOA([m2,m3,m4,m5],:).^2 - repmat(D_TDOA(m1,:).^2,4,1)];
                                [U,S,V]     = svd(Gamma_local);
                                Gamma_local = U(:,1:rank1)*S(1:rank1,1:rank1)*V(:,1:rank1)';
                                cnt([m1,m2,m3,m4,m5])= cnt([m1,m2,m3,m4,m5])+1;
                                KK(m1,cnt(m1),:) = Gamma_local(1,:);
                                KK(m2,cnt(m2),:) = Gamma_local(2,:);
                                KK(m3,cnt(m3),:) = Gamma_local(3,:);
                                KK(m4,cnt(m4),:) = Gamma_local(4,:);
                                KK(m5,cnt(m5),:) = Gamma_local(5,:);
                            end
                        end
                    end
                end
            end
            KK = KK(:,1:cnt(1),:);
            D_TDOA = squeeze(mean(KK,2));
        end;
        D_TDOA_new = D_TDOA;
    elseif dimension == 4  % 4-D case      
        rank1       = min([10 , N-1]);        
        for k = 1:iteration    
            cnt=zeros(M,1);
            KK = zeros(M,1000,N-1);
            for m1 = 1:M-5
                for m2 = m1+1:M-4
                    for m3 = m2+1:M-3
                        for m4 = m3+1:M-2
                            for m5 = m4+1:M-1
                                for m6 = m5+1:M
                                    Gamma_local = [D_TDOA([m1,m2,m3,m4,m5,m6],:) ; D_TDOA([m2,m3,m4,m5,m6],:).^2 - repmat(D_TDOA(m1,:).^2,5,1)];
                                    [U,S,V]     = svd(Gamma_local);
                                    Gamma_local = U(:,1:rank1)*S(1:rank1,1:rank1)*V(:,1:rank1)';
                                    cnt([m1,m2,m3,m4,m5,m6])= cnt([m1,m2,m3,m4,m5,m6])+1;
                                    KK(m1,cnt(m1),:) = Gamma_local(1,:);
                                    KK(m2,cnt(m2),:) = Gamma_local(2,:);
                                    KK(m3,cnt(m3),:) = Gamma_local(3,:);
                                    KK(m4,cnt(m4),:) = Gamma_local(4,:);
                                    KK(m5,cnt(m5),:) = Gamma_local(5,:);
                                    KK(m6,cnt(m6),:) = Gamma_local(6,:);
                                end
                            end
                        end
                    end
                end
            end
            KK = KK(:,1:cnt(1),:);
            D_TDOA = squeeze(mean(KK,2));
        end;            
        D_TDOA_new = D_TDOA;
    end
    
    %===== Step 5 of Algorithm: Refinement of TDOAs based on Theorem 2 ====
    Gamma_global = [D_TDOA_new ; D_TDOA_new(2:M,:).^2 - repmat(D_TDOA_new(1,:).^2,M-1,1)];
    [U,S,V]      = svd(Gamma_global);
    rank         = min([M+dimension , N-1]);
    Gamma_global = U(:,1:rank)*S(1:rank,1:rank)*V(:,1:rank)';
    D_TDOA_new   = Gamma_global(1:M,:);
    
    