Ok, I changed the algorithm a bit and now I think I nailed it and solved it for frequencies up to Nyquist and above! I'm just randomly choosing any sample in the period, from first to last and I don't have if there are consecutive samples, like the last sample from the previous period and first sample from the current one.
for(int j = 0; j < x->x_nchans; j++){
for(int i = 0, n = x->x_n; i < n; i++){
double hz = in1[j*n + i];
double step = hz * x->x_sr_rec; // phase step
step = step > 1 ? 1 : step < 0 ? 0 : step; // clip step
t_float imp = 0;
t_float r = x->x_rand[j]; // - step;
if(phase[j] >= r && ((lastphase[j] < r) || (x->x_1st[j])))
imp = velvet_random(x) > 0.5 ? 1 : -1;
out[j*n + i] = imp;
x->x_1st[j] = 0;
if(phase[j] >= 1.){
x->x_1st[j] = 1;
x->x_rand[j] = velvet_random(x);
while(phase[j] >= 1.)
phase[j] -= 1; // wrap phase
}
lastphase[j] = phase[j];
phase[j] += step;
}
}```