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[94],Mann[98] 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);

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.