I'm trying to modify some existing matlab code from 2 dimensions to 3 dimensions, and running in to a conceptual block. The original script generates fractals from 2d random noise by running an ifft(3) on the product of a power spectrum and random noise. In order to define the power spectrum I have to define the frequencies in U, V and W dimensions.
I'm having a mental block of what the W dimension frequencies should look like in relation to U and V. Obviously they wouldn't just be a transpose of U or V. Included a little extra code at the beginning and end if it helps clarify things. The var powerspectrum is currently returning a 2d matrix and it should be 3d. Thanks!
beta = -2;
imsize = [sidelength, sidelength, sidelength]; %length of sides of image in pixels;
% Generate the grid of frequencies. u is the set of frequencies along the first dimension
u = [(0:floor(imsize(1)/2)) -(ceil(imsize(1)/2)-1:-1:1)]'/imsize(1); % First quadrant are positive frequencies. Zero frequency is at u(1,1).
u = repmat(u,1,imsize(2)); % Reproduce these frequencies along ever row.
v = [(0:floor(imsize(2)/2)) -(ceil(imsize(2)/2)-1:-1:1)]/imsize(2); % v is the set of frequencies along the second dimension. For a square region it will be the transpose of u.
v = repmat(v,imsize(1),1); % Reproduce these frequencies along ever column.
w = [(0:floor(imsize(3)/2)) -(ceil(imsize(3)/2)-1:-1:1)]/imsize(3); % W is the set of frequencies along the *third* dimension.
w = repmat(w,imsize(3),1); % Reproduce these frequencies along every PAGE.
powerspectrum = (u.^2 + v.^2 + w.^2).^(beta/2); % Generate the power spectrum
phases = randn(imsize); % Generate a grid of random phase shifts
complexpattern = ifftn(powerspectrum.^0.5 .* (cos(2*pi*phases)+1i*sin(2*pi*phases))); % Inverse Fourier transform to obtain the the spatial pattern
Instead of using repmat you can reshape u, v and w into [sidelength 1 1] , [1 sidelength 1] and [1 1 sidelength] arrays respectively and let implicit expansion to do its work and create a [sidelength sidelength sidelength] array:
u = [(0:floor(imsize(1)/2)) -(ceil(imsize(1)/2)-1:-1:1)]'/imsize(1);
u = reshape(u,[],1)
v = [(0:floor(imsize(2)/2)) -(ceil(imsize(2)/2)-1:-1:1)]/imsize(2);
v = reshape(v,1,[]);
w = [(0:floor(imsize(3)/2)) -(ceil(imsize(3)/2)-1:-1:1)]/imsize(3);
w = reshape(w,1,1,[]);
powerspectrum = (u.^2 + v.^2 + w.^2).^(beta/2);
phases = randn(imsize);
complexpattern = ifftn(powerspectrum.^0.5 .* (cos(2*pi*phases)+1i*sin(2*pi*phases)));