function x_s_sigma = hierarchicalThresholdOperator(x, N, n, s, sigma)
X = reshape(x, [n, N]);
    
    % Blockwise threshold to keep sigma largest entries in each block
    X_sigma = zeros(size(X));   % This will hold the sigma-sparse blocks
    for blk = 1:N
        block = X(:, blk);
        % Find the indices of the largest "sigma" entries in magnitude
        [~, idx] = maxk(abs(block), sigma); 
        % Keep those sigma entries, zero out the rest
        blockSparse = zeros(1, n);
        blockSparse(idx) = block(idx);
        
        X_sigma(:, blk) = blockSparse;
    end
    
    % Compute the block norms (l2-norm) of each sigma-sparse block
    blockNorms = vecnorm(X_sigma);  % Size: N x 1
    
    % Keep the s blocks with the largest norms
    [~, blockIdx] = maxk(blockNorms, s);
    
    X_s_sigma = zeros(size(X_sigma));
    X_s_sigma(:, blockIdx) = X_sigma(:, blockIdx);
    
    % Reshape back to a length-(N*n) vector
    x_s_sigma = vec(X_s_sigma);
    % omega = find(abs(x_s_sigma) > 0); % for debugging purpose
end
