Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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:

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.