function h_hat = predict_channel_more_v_20(y, simPar, predict_Method, user,p)
    
    if predict_Method == "LS"
        h_hat = y;
%         h_hat = h_hat/norm(h_hat)*sqrt(simPar.N);
    elseif predict_Method == "MMSE_ISO"
%         mu = [0; 0; 10];% mu = [0.0001; -0.0001; 0.0001];
        mu = [0.0001; 0; 0.0001];
        theta = struct('rho', 2, 'v', [0; 0; 0], 'mu', mu);
        [Kha, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,theta, simPar.k0, simPar.c);
        [Khb, ~] = get_K(simPar.x,simPar.xf, simPar.p,simPar.pf, theta, simPar.k0, simPar.c);
%         eta = simPar.lambda/2/pi;
% n_azi = ( -(simPar.N-1)/2:1:(simPar.N-1)/2 )';
% % n_ele = ( -(N2-1)/2:1:(N2-1)/2 )';
% kernel_azi = exp(-0.1*pi*abs(simPar.BS.spacing*n_azi - simPar.BS.spacing*n_azi').^2);
% kernel_azi = besselj(0,0.85*abs(simPar.BS.spacing*n_azi - simPar.BS.spacing*n_azi')/eta);
% kernel_ele = exp(-0.096*pi*abs(d*n_ele - d*n_ele').^2/eta^2);
% Kernel_exp = kron(kernel_azi, kernel_azi);

%         h_hat =  kernel_azi*((kernel_azi + simPar.sigma2_n*eye(simPar.N))\y(33:64));
        h_hat =  Khb*((Kha + simPar.sigma2_n*eye(simPar.N*simPar.T))\y);

    elseif predict_Method == "GPR"
        if simPar.channelEstimator.GPR.mixedKernel == false
            initialTheta = struct( 'v',    [0;0;0], ...
                                'mu',   [0.001; 0; 0.001]);%[-0.2;+0.000; -0.1]
        
            [~, optTheta] = fit_GPR_v_20(y, initialTheta, simPar, false);
           [optKh, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,optTheta, simPar.k0, simPar.c);
            [optKhb, ~] = get_K(simPar.x,simPar.xf, simPar.p,simPar.pf, optTheta, simPar.k0, simPar.c);

        else
            initialTheta = struct( 'v',    [rand();0;rand()], ...
                                'mu',   [0.001; 0; 0.001]);
            [~, optTheta] = fit_GPR_v_20(y, initialTheta, simPar, false);
            optKh = zeros(simPar.N*simPar.T);
            optKhb = zeros(simPar.N*simPar.T_pre,simPar.N*simPar.T);
            for idx = 1:simPar.channelEstimator.GPR.mixKernelNum
                [K, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,optTheta.kernelParams(idx), simPar.k0, simPar.c);
                [Kb, ~] = get_K(simPar.x,simPar.xf, simPar.p,simPar.pf, optTheta.kernelParams(idx), simPar.k0, simPar.c);
                optKh = optKh + optTheta.kernelWeights(idx) * K;
                optKhb = optKhb+optTheta.kernelWeights(idx) * Kb;
            end
        end
        h_hat = optKhb*((optKh + simPar.sigma2_n*eye(simPar.N*simPar.T))\y); 
    
    elseif predict_Method == "GPRmu"
        mu = 1*user;
        theta = struct('rho', 2, 'v', [0; 0; simPar.v], 'mu', mu);
        [Kha, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,theta, simPar.k0, simPar.c);
        [Khb, ~] = get_K(simPar.x,simPar.xf, simPar.p,simPar.pf, theta, simPar.k0, simPar.c);
        h_hat =  Khb*((Kha + simPar.sigma2_n*eye(simPar.N*simPar.T))\y);

%         if simPar.channelEstimator.GPR.mixedKernel == false
%             initialTheta = struct( 'v',    [0;0;20], ...
%                                 'mu',   [20; 0; 5]);%[-0.2;+0.000; -0.1]
%         
%             [~, optTheta] = fit_GPR(y, initialTheta, simPar, false);
%             [optKh, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,optTheta, simPar.k0, simPar.c);
%             [optKhb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, optTheta, simPar.k0, simPar.c);
% 
%         else
%             initialTheta = struct( 'v',    [0;0;1], ...
%                                 'mu',   user);
%             [~, optTheta] = fit_GPR(y, initialTheta, simPar, false);
%             optKh = zeros(simPar.N);
%             for idx = 1:simPar.channelEstimator.GPR.mixKernelNum
%                 [K, ~] = get_K(simPar.x,simPar.x, simPar.p, optTheta.kernelParams(idx), simPar.k0, simPar.c);
%                 optKh = optKh + optTheta.kernelWeights(idx) * K;
%             end
%         end
%         h_hat = optKhb*((optKh + simPar.sigma2_n*eye(simPar.N*simPar.T))\y);

    elseif predict_Method == "AMP"
        A = dftmtx(simPar.N)/sqrt(simPar.N); 
        
        x_hat = camp2(y/sqrt(simPar.N), A, simPar.channelEstimator.AMP.lambda);
        h_hat = A*x_hat*sqrt(simPar.N); 
    
    elseif predict_Method == "OMP"
        A = dftmtx(simPar.N); 
        
        x_hat = OMP(y, A, 7);
        h_hat = A*x_hat; 

    elseif predict_Method == "GPRbes"
        [h_hat, ~] = fit_GPRbes(y, simPar); 
    
    elseif predict_Method == "Prony"
        R = y(:,1:simPar.T-1)'*y(:,1:simPar.T-1)/(simPar.N);
        v = y(:,1:simPar.T-1)'*y(:,simPar.T)/(simPar.N);
        P2 = inv(R)*v;
%         p = pinv(y(:,1:simPar.T-1))*y(:,simPar.T);
        h_hat = y(:,2:simPar.T)*p;
% h_hat = h_hat/norm(h_hat)*simPar.N;
    elseif predict_Method == "AR"
        fmax = simPar.v/simPar.lambda;
        T = simPar.Ts;
        z = 2*pi*fmax*T;
%         %一阶
%         aa = besselj(0, z);
%         h_hat = aa*y;
%         %二阶
%         Rxx=[1,besselj(0, z);besselj(0, z),1];
%         v_ = [besselj(0, z);besselj(0, 2*z)];
%         a_hat = inv(Rxx)*v_;
%         h_hat = a_hat(1)*y(1:simPar.N)+a_hat(2)*y(simPar.N+1:2*simPar.N);
        %T阶
        Rxx = toeplitz(besselj(0, (0:simPar.T-1)*z));
        Rxx = toeplitz(exp(-(0:simPar.T-1)*z));
        mu = [0.0001; 0; 0.0001];
        theta = struct('rho', 2, 'v', [0; 0; 20], 'mu', mu);
        [Kha, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,theta, simPar.k0, simPar.c);
        [Khb, ~] = get_K(simPar.x,simPar.xf, simPar.p,simPar.pf, theta, simPar.k0, simPar.c);
%         Rxx = zeros(simPar.T,simPar.T,'double');
        
        h_hat = zeros(simPar.N, 1);
        for jd=1:simPar.N
        for id=1:simPar.T
            r = zeros(1,simPar.T,'double');
            v_ = zeros(1,simPar.T,'double');
            r(id) = Kha(jd,(id-1)*simPar.N+jd);
            v_(id) = Khb(jd,(id-1)*simPar.N+jd);
        end
%         Rxx = toeplitz(r);
        y2 = [y(1:simPar.N),y(simPar.N+1:simPar.N*2)];
%         Rxx = y(1:simPar.N)'*y(1:simPar.N);

        v_ = [besselj(0, (1:simPar.T)*z)];
        v_ = [exp(-(1:simPar.T)*z)];
        v_ = v_.';
%         v_ = y(1:simPar.N)'*y(simPar.N+1:simPar.N*2);

%         a_hat = Rxx\v_;
        a_hat = pinv(Rxx)*v_;
        
        mx_number = simPar.N;
% mz_number = 50;
% my_number = 50;
% threshold = 0.2;
% threhold = 1e-3;
% mulength = 40;
% mz_total = -mz_number:mz_number;
% my_total = -my_number:my_number;
lam = simPar.lambda;
dy = lam/8; dz =  dy;
% iter_test = 30;
% R = zeros(2*my_number+1,2*mz_number+1);
% N = (2*mz_number+1)*(2*my_number+1);

% bsPosition = [22.287495, 114.140706]; % lat, lon
% bsAntSize = [mx_number 1];                    % number of rows and columns in rectangular array (base station)
% bsArrayOrientation = [0 0].';       % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards) in deg
% uePosition = [22.287551+rand(1)*0.000008,114.140702+rand(1)*8*0.000001];  % lat, lon
% ueAntSize = [1 1];                    % number of rows and columns in rectangular array (UE).
% ueArrayOrientation = [180 0].';      % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards)  in deg
% reflectionsOrder = 2;                 % number of reflections for ray tracing analysis (0 for LOS)

bsPosition = [22.287+rand(1)*0.001, 114.140+rand(1)*0.001]; % lat, lon
bsAntSize = [mx_number 1];                    % number of rows and columns in rectangular array (base station)
bsArrayOrientation = [90 0].';       % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards) in deg
uePosition = [bsPosition(1)-0.000001*4+rand(1)*0.000008,bsPosition(2)+0.000001*20+rand(1)*8*0.000001-0.000001*4];  % lat, lon
ueAntSize = [1 1];                    % number of rows and columns in rectangular array (UE).
ueArrayOrientation = [270 0].';      % azimuth (0 deg is East, 90 deg is North) and elevation (positive points upwards)  in deg
reflectionsOrder = 2; 

bsSite = txsite("Name","Base station", ...
    "Latitude",bsPosition(1),"Longitude",bsPosition(2),...
    "AntennaAngle",bsArrayOrientation(1:2),...
    "AntennaHeight",4,...  % in m
    "TransmitterFrequency",3e8/lam);

ueSite = rxsite("Name","UE", ...
    "Latitude",uePosition(1),"Longitude",uePosition(2),...
    "AntennaHeight",4,... % in m
    "AntennaAngle",ueArrayOrientation(1:2));

% bsSite = txsite("Name","Base station", ...
%     "Latitude",bsPosition(1),"Longitude",bsPosition(2),...
%     "AntennaHeight",0,...  % in m
%     "TransmitterFrequency",3e8/lam);
% 
% ueSite = rxsite("Name","UE", ...
%     "Latitude",uePosition(1),"Longitude",uePosition(2),...
%     "AntennaHeight",0);

% show(bsSite);
% show(ueSite);

pm = propagationModel("raytracing","Method","sbr","MaxNumReflections",reflectionsOrder);
rays = raytrace(bsSite,ueSite,pm,"Type","pathloss");

% plot(rays{1})

pathToAs = [rays{1}.PropagationDelay]-min([rays{1}.PropagationDelay]);  % Time of arrival of each ray (normalized to 0 sec)
avgPathGains  = -[rays{1}.PathLoss];                                    % Average path gains of each ray
pathAoDs = [rays{1}.AngleOfDeparture];                                  % AoD of each ray
pathAoAs = [rays{1}.AngleOfArrival];                                    % AoA of each ray
isLOS = any([rays{1}.LineOfSight]);                    
        cdl_model = nrCDLChannel('ChannelFiltering', false);
cdl_model.DelayProfile = 'Custom';
cdl_model.PathDelays = pathToAs;
cdl_model.AveragePathGains = avgPathGains;%!
cdl_model.AnglesAoD = pathAoDs(1,:);       % azimuth of departure
cdl_model.AnglesZoD = 90-pathAoDs(2,:);    % channel uses zenith angle, rays use elevation
cdl_model.AnglesAoA = pathAoAs(1,:);       % azimuth of arrival
cdl_model.AnglesZoA = 90-pathAoAs(2,:);    % channel uses zenith angle, rays use elevation
cdl_model.HasLOSCluster = isLOS;
simPar.sta = 350;
%     R = zeros(simPar.T,simPar.T);
%     DD = zeros(simPar.T+1,simPar.sta);
%     v = zeros(simPar.T,1);
%     simPar.a = zeros(simPar.N,simPar.T);
%     D = zeros(simPar.N,simPar.T+1,simPar.sta);
%     Num = zeros(simPar.T+1,simPar.sta);
        cdl_model = nrCDLChannel('ChannelFiltering', false);
        cdl_model.MaximumDopplerShift       = simPar.v/simPar.lambda;
        cdl_model.UTDirectionOfTravel       = [0;0];
        cdl_model.CarrierFrequency          = simPar.fc;
        cdl_model.SampleRate                = 1/simPar.Ts;
        cdl_model.SampleDensity             = Inf;
        cdl_model.NumTimeSamples            = simPar.T+simPar.sta;
        cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
        cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
        cdl_model.Seed                      = randi(10000);
        [pathGains, ~] = cdl_model();
        for idx1 = 1:simPar.T+simPar.sta
             h1(:,idx1) = reshape(sum(pathGains(idx1, :, :, :), 2), [simPar.N, 1]);
        end
        for j = 1:simPar.T+1
        simPar.h((j-1)*simPar.N+1:j*simPar.N) = h1(simPar.sta+j-1);
        end
        for idx2 = 1:simPar.T+1
            for idx3 = 1:simPar.sta
           D(:,idx2,idx3) = h1(:,idx3+idx2-1);
           Num(idx2,idx3) = idx3+idx2-1;
            end
        end
    for idx = 1:simPar.N
    DD(:,:) = D(idx,:,:);
    RQ = DD*DD'/simPar.sta;
    R = RQ(1:simPar.T,1:simPar.T)+eye(simPar.T)*0.01;
    v = flip(RQ((simPar.T+1),(1:(simPar.T))));
    v = v.';
%     [V,D2,W] = eig(R)
    simPar.a(idx,:) = inv(R)*v;
    end
        for i = 1:simPar.T
%             h_hat = h_hat+a_hat(i)*y((simPar.T-i)*simPar.N+1:(simPar.T-i+1)*simPar.N);
%             h_hat(jd) = h_hat(jd)+a_hat(i)*y((simPar.T-i)*simPar.N+jd);
            h_hat(jd) = h_hat(jd)+simPar.a(jd,i)*y((simPar.T-i)*simPar.N+jd);
        end

        end
%  v_ = [exp(-(1:simPar.T)*z)];
%         v_ = v_.';
% 
% %         a_hat = Rxx\v_;
%         a_hat = inv(Rxx)*v_;
%         h_hat = zeros(simPar.N, 1);
% 
%         for i = 1:simPar.T
%             h_hat = h_hat+a_hat(i)*y((simPar.T-i)*simPar.N+1:(simPar.T-i+1)*simPar.N);
%         end
%         a_hat = zeros(simPar.N,2);
%         h_hat = zeros(simPar.N,1);
%         p = 2;%order
%         for i=1:simPar.N
%             Rxx = zeros(p, p); % 自相关矩阵
%             v = zeros(p, 1); % Yule-Walker方程的向量
%             % 计算自相关矩阵和Yule-Walker方程的向量
%             for j = 1:p
%                 for k = 1:p
%                 if k == j
%                     Rxx(lag+1, lag+1) = abs(y(i)).^2;
%                 else
%                     m = 0;
%                     for ii = 1:p
%                         m = m+abs(y(i+(k-1)*simPar.N)).*conj(H(:, lag+1:end));
%                     end   
%                     Rxx(j, k) = sum(abs(y(i+(k-1)*simPar.N)).*conj(H(:, lag+1:end)))/p;
%                 end
%                 end
%             end
%                 for k = 1:p
%                     if lag == 0
%                         v(k) = sum(abs(H(:, k)).*conj(H(:, 1)));
%                     else
%                         v(k) = sum(sum(abs(H(:, k+lag:end)).*conj(H(:, 1:end-lag))));
%                     end
%                 end
%             
%         end

%         for i=1:simPar.N
%             [a_hat(i,:)] = ar([y(i),y(simPar.N+i)],2,'ls');
%             h_hat(i) = a_hat(i,1)*y(i)+a_hat(i,2)*y(simPar.N+i);
%         end

    else
        error("Error: Channel estimation method %s is not implemented.", predict_Method);
    end

end

