function [NMSE, pilotOverhead] = estimateChannel(H, Correlation, MIMO_params, method)
    
    sigma2 = MIMO_params.sigma2; 
    NTx = MIMO_params.TxArray.N; 
    NRx = MIMO_params.RxArray.N; 
    Npilot = MIMO_params.Npilot; 

    if method == "OMP"
        


    elseif method == "PSWF"
        % Step 1: Support detection

        Hk = (MIMO_params.RxArray.dict)'*H*(MIMO_params.TxArray.dict); 
        
        figure; 
        imagesc(abs(Hk)); axis equal; colorbar; 
        

        

        fprintf('fuck\n'); 

    elseif method == "PSWFPilotMMSE"
        X = pilotGeneration(MIMO_params, "PSWF"); 
        N = (randn([NRx, Npilot]) + 1j*randn([NRx, Npilot]))/sqrt(2); 
        
        targetSNR = MIMO_params.P_channelEst;   % SNR observed on each receiving antenna. 
        sigma2 = (norm(H, 'fro')^2/NRx)/targetSNR; 
        eSNR = norm(H, 'fro')^2/(NRx*NTx) / sigma2;

        Y = H*X+sqrt(sigma2)*N; 

        Hhat = zeros(NRx, NTx); 
        % Model: y = Ax + w
        for idx = 1:NRx
            Hhat(idx, :) = ((X.')'*(X.')+eye(NTx)/eSNR)\((X.')'*(Y(idx, :)).'); 
        end

    elseif method == "RandomPilotMMSE"
        X = pilotGeneration(MIMO_params, "Random"); 
        N = (randn([NRx, Npilot]) + 1j*randn([NRx, Npilot]))/sqrt(2); 
        
        
        % control receive SNR. 
        targetSNR = MIMO_params.P_channelEst;   % SNR observed on each receiving antenna. 
        sigma2 = (norm(H, 'fro')^2/NRx)/targetSNR; 
        eSNR = norm(H, 'fro')^2/(NRx*NTx) / sigma2;

        Y = H*X+sqrt(sigma2)*N; 

        Hhat = zeros(NRx, NTx); 
        % Model: y = Ax + w
        for idx = 1:NRx
            Hhat(idx, :) = ((X.')'*(X.')+eye(NTx)/eSNR)\((X.')'*(Y(idx, :)).'); 
        end
        
    
    elseif method == "RandomPilotMMSEwithPrior"
        X = pilotGeneration(MIMO_params, "Random"); 
        N = sqrt(sigma2)*(randn([NRx, Npilot]) + 1j*randn([NRx, Npilot]))/sqrt(2); 
        Y = H*X+N; 

        % Model: H = (Rr)^{1/2} * W * (Rt)^{1/2,T} with W standard complex
        % Gaussian. 
        
        halfRr = Correlation.halfRr; 
        halfRt = Correlation.halfRt; 
        Dr = MIMO_params.RxArray.dict;
        Dt = MIMO_params.TxArray.dict; 

        A = (2/MIMO_params.K)*(Dr)*halfRr; 
        B = (2/MIMO_params.K)*halfRt.'*(Dt')*X; 

        % temp = 1/sigma2*(A'*Y*B'); 
        % what = (eye(NTx*NRx)+(1/sigma2)*kron((B*B').', (A'*A)))\temp(:); 
        % What = reshape(what, [NRx, NTx]); 
        
        What = ((sigma2*eye(MIMO_params.K)+A'*A)\(A'*Y*B'));
        What = What/(sigma2*eye(MIMO_params.K) + B*B'); 

        Hhat = Dr * halfRr * What * (halfRt).' * Dt' *(2/MIMO_params.K)^2; 
    
    else
        error('Method %s not implemented.', method); 

    end

    NMSE = norm(Hhat - H, 'fro')^2/(norm(H, 'fro')^2); 
    pilotOverhead = Npilot; 

end

