# Keep the exact frequency content - Multidimensional fft

In Matlab I've got to generate a Gaussian process embodying a certain frequency content. Not to go multidimensional, I have used reshape function. I explain better with a piece of code:

dk1 = 2*pi*L/L1;
dk2 = 2*pi*L/L2;
dk3 = 2*pi*L/L3;
k1 = dk1*(-N1/2:(N1/2-1));
k2 = dk2*(-N2/2:(N2/2-1));
k2min = min(k2);
k2max = max(k2);
k3 = dk3*(-N3/2:(N3/2-1));
k3min = min(k3);
k3max = max(k3);
k = allcomb(k1,k2,k3);
k1 = reshape(k(:,1),[1,N]);
k2 = reshape(k(:,2),[1,N]);
k3 = reshape(k(:,3),[1,N]);
k = sqrt(k1.^2 + k2.^2 + k3.^2);
idx = (k1~=0);


allcomb may be found here: http://www.mathworks.com/matlabcentral/fileexchange/10064

Afterwards, I proceed as follows:

k30 = k3(idx) + beta.*k1(idx);
k0 = sqrt(k1(idx).^2 + k2(idx).^2 + k30.^2);

C1 = (beta.*k1(idx).^2.*(k1(idx).^2 + k2(idx).^2 - k3(idx).*k30))./(k(idx).^2.*(k1(idx).^2 + k2(idx).^2));
C2 = ((k2(idx).*k0.^2)./((k1(idx).^2 + k2(idx).^2).^(3/2))).*atan2((beta.*k1(idx).*sqrt(k1(idx).^2 + k2(idx).^2)),(k0.^2 - k30.*k1(idx).*beta));

xhsi1 = C1 - (k2(idx)./k1(idx)).*C2;
xhsi2 = (k2(idx)./k1(idx)).*C1 + C2;

Ek0 = 1.453.*k0.^4./((1 + k0.^2).^(17/6));

B = sigma_iso*sqrt((2*pi.^2*L.^3.*Ek0)./(V.*k0.^4));

%% C-matrix Calculation for all those k satisfying the condition |k| > 3

C = zeros(3,3,N);
% In this case, the approximation proposed in Mann,Mann and IEC
% 61400-1 3rd Ed. will be used
C(1,1,idx) = B.*(k2(idx).*xhsi1);
C(1,2,idx) = B.*(k3(idx) - k1(idx).*xhsi1 + beta.*k1(idx));
C(1,3,idx) = B.*(-k2(idx));

C(2,1,idx) = B.*(k2(idx).*xhsi2 - k3(idx) - beta.*k1(idx));
C(2,2,idx) = B.*(-k1(idx).*xhsi2);
C(2,3,idx) = B.*(k1(idx));

C(3,1,idx) = B.*(k2(idx).*((k0./k(idx)).^2));
C(3,2,idx) = B.*(-k1(idx).*((k0./k(idx)).^2));

nn = normrnd(0,1,[3,N]);

u = zeros(3,N);
for l = 1:size(C,3)
u(:,l) = C(:,:,l)*nn(:,l);
end


Then, I need a 3D fft of u; thus, I use reshape once more as

u = reshape(u,[3,N1,N2,N3]);
U = fftn(u);


But I believe that this flow is somehow affecting and corrupting the right frequency content.

Generally, I would go for:

[k1,k2,k3] = ndgrid(k1,k2,k3);
k = sqrt(k1.^2 + k2.^2 + k3.^2);


further to u calculation. The core of the code would switch to 5D instead of 3D.

Do you have any idea how to improve this code? Or may you spot the flaw?

Any support or criticism are welcome.