function [W, V] = getMIMOCombinerPrecoder_new(MIMO_params, IntervalsRx, IntervalsTx, Ncombiner, method, thres)
    % W: Combiner matrix; V: Precoder matrix. 

    NRx = MIMO_params.RxArray.N; 
    NTx = MIMO_params.TxArray.N; 

    if method == "PSWF"
        cntBasisRx = 0; 
        % selDSSRx = zeros(NRx, 1); 
        
        LambdasRx = zeros(NRx*length(IntervalsRx), 1); 
        DSSRx_vectors = zeros(NRx, NRx*length(IntervalsRx)); 
        
        % Gather all the DPSSs 
        for idx = 1:length(IntervalsRx)
            kx = IntervalsRx{idx}(1); 
            rx = IntervalsRx{idx}(3); 
    
            [DSSRx, LambdaRx] = getDiscreteSlepianSequence(NRx, kx-rx, kx+rx); 
            
            DSSRx_vectors(:, (idx-1)*NRx+1:idx*NRx) = DSSRx; 
            LambdasRx((idx-1)*NRx+1:idx*NRx) = LambdaRx; 

            % for ii = 1:NRx
            %     if LambdaRx(ii) >= thres
            %         selDSSRx(:, cntBasisRx+1) = DSSRx(:,ii); 
            %         cntBasisRx = cntBasisRx + 1;
            %         LambdasRx = [LambdasRx, LambdaRx(ii)]; 
            %     end
            % end
        end
        
        % Select the top Ncombiner DPSSs 
        [~, order] = sort(LambdasRx, "descend"); 
        selDSSRx = zeros(NRx, Ncombiner);  

        if Ncombiner < NRx*1
        selDSSRx(order(1:Ncombiner)) = true; 
        end



        cntBasisTx = 0; 
        selDSSTx = zeros(NTx, 1); 
        
        LambdasTx = []; 
        for idx = 1:length(IntervalsTx)
            kx = IntervalsTx{idx}(1); 
            rx = IntervalsTx{idx}(3); 
    
            [DSSTx, LambdaTx] = getDiscreteSlepianSequence(NTx, kx-rx, kx+rx); 
            
            for ii = 1:NTx
                if LambdaTx(ii) >= thres
                    selDSSTx(:, cntBasisTx+1) = DSSTx(:,ii); 
                    cntBasisTx = cntBasisTx + 1;
                    LambdasTx = [LambdasTx, LambdaTx(ii)]; 
                end
            end
        end

        if cntBasisRx >= Ncombiner
            [~, order] = sort(LambdasRx, 'descend'); 

            p = randperm(cntBasisRx, Ncombiner); 
            % p = order(1:Ncombiner); 
            W = selDSSRx(:, p); 
        else
            % randX = (randn([NRx, Ncombiner-cntBasisRx]) + 1j*randn([NRx, Ncombiner-cntBasisRx]))/sqrt(2); 
            % randX = randX ./ vecnorm(randX); 
            randChosen = randi([1, cntBasisRx], [1, Ncombiner-cntBasisRx]); 
            W = [selDSSRx, selDSSRx(:, randChosen)]; 
            % W = W(:, randperm(Ncombiner)); 
        end
    
        if cntBasisTx >= Ncombiner
            [~, order] = sort(LambdasRx, 'descend'); 

            p = randperm(cntBasisTx, Ncombiner); 
            % p = order(1:Ncombiner); 
            V = selDSSTx(:, p);

        else
            randChosen = randi([1, cntBasisTx], [1, Ncombiner-cntBasisTx]); 
            V = [selDSSRx, selDSSTx(:, randChosen)]; 
            % W = W(:, randperm(Ncombiner)); 
        end

    elseif method == "Random"
        randX = (randn([NRx, Ncombiner]) + 1j*randn([NRx, Ncombiner]))/sqrt(2); 
        W = randX ./ vecnorm(randX); 

        randY = (randn([NTx, Ncombiner]) + 1j*randn([NTx, Ncombiner]))/sqrt(2); 
        V = randY ./ vecnorm(randY); 
        
    else
        error('Method %s does not exist. ', method); 
    end

end
