How do I normalize the signal generated from wfbm in matlab?

Question

I am using matlab's wavelet fractional Brownian motion function in order to generate 1D point-like data of a diffusive particle in the physical regimes: sub-diffusion, super-diffusion and normal diffusion.

The problem I encounter with is that the time normalization/variance is weird.

For example for Hurst parameter equals 0.5 (regular Brownian motion) I get standard deviation which isn't unity (1):

>> std(diff(wfbm(0.5,1e6)))

ans =

    0.3955

Due to the above, I am not sure how to re-normalize all the 3 trajectories I create for the 3 diffusion cases (sub, super, normal).

---

I generated trajectories for N pointlike particles of length M:

M=500;
N=200;
nd  = zeros(M,N);
sub = zeros(M,N);
sup = zeros(M,N);
Hsub = 0.25;
Hsup = 0.75;

for j=1:N
    nd(:,j)  = wfbm(0.5, M, 15, 'db10');
    sub(:,j) = wfbm(Hsub,M, 10, 'db10');
    sup(:,j) = wfbm(Hsup,M, 10, 'db10');
end

---

Here is how function is implemented in matlab and generates the signal, however I am not sure how to modify it to have a proper brownian motion:

tmp = conv(randn(1,len+nbmax),ckbeta);
tmp = cumsum(tmp);
CA  = wkeep(tmp,len,'c');
for j=0:nblev-1
    CD  = 2^(j/2)*4^(-s)*2^(-j*s)*randn(1,len);
    len = 2*len-nbmax;
    CA  = idwt(CA,CD,fs1,gs1,len);
end
fBm  = wkeep(CA,L,'c');
fBm  = fBm-fBm(1);

I was trying to understand it from the paper which says it's possible to control the variance of fBm:

This is citation 7 from the snapshot above.