function [ est_H,est_mat_MS,est_mat_BS,est_mat_F,est_gain,est_theta_1,est_theta_2,est_theta_3 ] = Par_estimation(Y,My,Mz,A,B,C,V,F,S,dim_ratio,zl)
     
     M = size(V,1);
     Nt = size(F,1);
     P0 = size(S,2);

     [est_theta_1] = extraction_f_UPA(My,Mz,V.',A,60);
     [est_theta_2] = extraction_f_ULA(Nt,F.',B,dim_ratio);
     est_theta_3 = -log(zl)/1j;
     
     est_mat_MS=steering_matrix_UPA(My,Mz,est_theta_1);
     est_mat_BS=steering_matrix_ULA(Nt,est_theta_2);
     est_mat_F=steering_matrix_tao(P0,est_theta_3);

     Equal_A=khatrirao(F.'*est_mat_BS,V.'*est_mat_MS); 
     Equal_y=tenmat(Y,3);
     Equal_y=conj(Equal_y.data);
     e_C = Equal_y*pinv(Equal_A.');
     est_gain = scale(e_C,S*est_mat_F);    
     
     est_H=ktensor(est_gain,est_mat_MS,est_mat_BS,S*est_mat_F);
     est_H=tensor(est_H);
     
      function f = extraction_f_UPA(My,Mz,R,Ob,dictionary_granularity)
            nn1 = My*dictionary_granularity;
            nn2 = Mz*dictionary_granularity;
            grid1 = ((1:nn1)*2*pi)/nn1;
            grid2 = ((1:nn2)*2*pi)/nn2;
            d1 = exp(1j*(0:(My-1))'*grid1);
            d2 = exp(1j*(0:(Mz-1))'*grid2);
            D = kron(d1,d2);
            D = R*D;
            [~,index] = max(abs(normalization_column(D)'*normalization_column(Ob)));
            f2 = mod(index,nn2);
            f1 = (index-f2)/nn2+1;
            if ~isempty(find(f2==0, 1))
                [~,indx] = find(f2==0);
                f2(indx)=nn2;
                f1(indx)=f1(indx)-1;
            end
            f = [grid1(f1).',grid2(f2).'];
      end

    function f = extraction_f_ULA(N,R,Ob,dictionary_granularity)
            nn = N*dictionary_granularity;
            grid = ((1:nn)*2*pi)/nn;
            D = exp(1j*(0:(N-1))'*grid);
            D = R*D;
            [~,index] = max(abs(normalization_column(D)'*normalization_column(Ob)));
            f = grid(index).';
    end

    function [A,lambda] = normalization_column(A) 
        lambda = sqrt(sum(abs(A).^2,1));
        M = size(A,1);
        A = A./(ones(M,1)*lambda);
    end

    function lambda=scale(A,B)
        D=pinv(B)*A;
        lambda=diag(D);
    end
end

