So far I haven't found a way to compute the circular convolution (of size N for sequences greater than N) that doesn't involve applying linear convolution + aliasing. Which is my original question

However, I have found in Sophocles J. Orfanidis' Introduction to signal processing page 518 a clever algorithm that allows me to obtain circular convolution from linear convolution + aliasing in an understandable way.

Here's my explanation:

First we must obtain linear convolution, then we place linear convolution output on time axis at intervals of L samples (sort of a periodic extension). As linear convolution output is greater than L samples, some samples will overlap (actually samples from more than one "period" will overlap on the main "period" from 0 to L-1). In obtaining circular convolution we are only interested in the interval 0 to L-1 (or 1 to L in Matlab)

For instance this example from the book computing 3 points circular convolution:

[1, 3, 3, 5, 3, 7, 4, 3, 3, 0, 1]  main period (direct linear convolution output)
[5, 3, 7, 4, 3, 3, 0, 1, 0, 0, 0]  1 period to the left contribution on main period
[4, 3, 3, 0, 1, 0, 0, 0, 0, 0, 0]  2 periods to the left contribution on main period
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]  3 periods to the left contribution on main period

No further periods to the left are considered as they don't generate aliasing in the main period (linear convolution output is smaller than 3 shifts of 3 samples). In the end to form the sequence the first 3 positions (of main period) are summed in order to form the 3 point circular convolution which is (1+5+4+0,3+3+3+1,3+7+3+0) = (10,10,13)

The periods to the right are ignored as they don't contribute to aliasing in the main "period"

Here's my Matlab code (for generic sequences):

modulo = 4;
sequence1  = [0.5,1,1,0.5];
sequence2  = [0.5,1,1,0.5];
linealConvolution = conv(sequence1,sequence2);

stem(linealConvolution)
hold on

periods = linealConvolution;
times = 1;

while length(linealConvolution) >= times*modulo
        
    %shift the linealConvolution sequence modulo(amount) samples at a time and add
    %it to the previosly stored sequence (this adds the contribution of time 
    %aliased lineal convolutions on adjacent periods)
    
    aliasingPeriodToTheLeft = [ linealConvolution(times*modulo + 1: end) zeros(1,times*modulo) ] ;
    stem(aliasingPeriodToTheLeft)
    
    periods = periods + aliasingPeriodToTheLeft;
    times = times+1;
    
end

inTheEnd = periods(1:modulo);


figure
stem(inTheEnd)

%checking
isequal(cconv(sequence1,sequence2,modulo),inTheEnd)