Random (in-place) shuffle of input wave's elements

function shuffle(inwave)    //  in-place random permutation of input wave elements
    wave inwave
    variable N  =   numpnts(inwave)
    variable i, j, emax, temp
    for(i = N; i>1; i-=1)
        emax = i / 2
        j =  floor(emax + enoise(emax))     //  random index
//      emax + enoise(emax) ranges in random value from 0 to 2*emax = i
        temp        = inwave[j]
        inwave[j]       = inwave[i-1]
        inwave[i-1] = temp
    endfor
end
 
function test(npts, nshuffles)
    variable npts, nshuffles
    make/O/N=(nshuffles, npts) wave0 = q+1  //  2D  (trials) x (array size)
    make/O/N=(npts) wtemp                   //  temporary column (trial) wave
    variable i
    for(i=0; i<nshuffles; i+=1)                 //  shuffle each row
        wtemp = wave0[i][p]
        shuffle(wtemp)
        wave0[i][] = wtemp[q]
    endfor
    Make/O/N=(npts) wmean   //  mean of the trials column for each array 'bin'
    for(i=0; i<npts; i+=1)
        MatrixOp/O wdest = col(wave0,i)
        wmean[i] = mean(wdest)  //  in this test, should approach (npts+1)/2
    endfor
end