How can I spped up my shuffle test for a large array?

I'm working on shuffling a large array(Jy) which is about 400 (numtrial)x 5(numtaper) x 80(numf) x 96(npair). I need to shuffle the first dimension for each column. My way is: 1.generate the index array by randperm, e.g. index is a 400 x 500(nPerm) matrix; 2.repmat the index array up to 400 x 5 x 80 x 96 x 500; 3.repmat the original array up to 400 x 5 x 80 x 96 and index the original large array by the index array to fully shuffle it. But it costs a large amount of time at the third step which is to index a huge array. Do you have any idea to improve this and save time?

%first step
I = zeros(numtrial,1,1,1,nPerm);
for iP = 1:nPerm
    I(:,1,1,1,iP) = randperm(numtrial);
end
%second step
i = numtrial*((1:1:(numtaper*numf*npair*nPerm))-1);
I = repmat(I,1,numtaper,numf,npair,1)+reshape(i,1,numtaper,numf,npair,nPerm);
%third step
Jyrand = repmat(Jy,1,1,1,1,nPerm);
Jyrand = Jyrand(I);%this line costs too much time