Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

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.

share|improve this question

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.