Directional artifacts in MATLAB randn arrays?

Question

I'm generating 3d fractal noise in MATLAB using a variety of methods. It's working relatively well, but I'm having an issue where I see vertical striping artifacts in my noise. This happens regardless of what data type or resolution I use.

Edit: I figured it out. The solution is posted as an answer below. Thanks everyone for your thoughts and guidance!

expo = 2^6;
dims = [expo,expo,expo];
beta = -4.5;

render = randnd(beta, dims); % Create volumetric fractal
render = render - min(render); % Set floor to zero
render = render ./ max(render); % Set ceiling to one
%render = imbinarize(render); % BW Threshold option
render = render .* 255; % For greyscale

slicer = 1; % Turn on image slicer/saver
i = 0; % Page counter
format = '.png';
imagename = '___testDump/slice';

imshow(render(:,:,1),[0 255]); %Single test image

if slicer == 1
    for c = 1:length(render)
    i = i+1;
    pagenumber = num2str(i);
    filename = [imagename, pagenumber, format];
    imwrite(uint8(render(:,:,i)),filename)
    end
end

function X = randnd(beta,varargin)

seed = 999;
rng(seed); % Set seed

%% X = randnd(beta,varargin)

% Based on similar functions by Jon Yearsley and Hristo Zhivomirov
% Written by Marcin Konowalczyk
% Timmel Group @ Oxford University
 
%% Parse the input
narginchk(0,Inf); nargoutchk(0,1);
 
if nargin < 2 || isempty(beta); beta = 0; end % Default to white noise
assert(isnumeric(beta) && isequal(size(beta),[1 1]),'''beta'' must be a number');
assert(-6 <= beta && beta <= 6,'''beta'' out of range'); % Put on reasonable bounds
 
%% Generate N-dimensional white noise with 'randn'
X = randn(varargin{:});
if isempty(X); return; end; % Usually happens when size vector contains zeros
 
% Squeeze prevents an error if X has more than one leading singleton dimension
% This is a slight deviation from the pure functionality of 'randn'
X = squeeze(X);
 
% Return if white noise is requested
if beta == 0; return; end;
 
%% Generate corresponding N-dimensional matrix of multipliers
N = size(X);
% Create matrix of multipliers (M) of X in the frequency domain
M = []; 
for j = 1:length(N)
    n = N(j);
    
    if (rem(n,2)~=0) % if n is odd
        % Nyquist frequency bin does not show up in odd-numbered fft
        k = ifftshift(-(n-1)/2:(n-1)/2);
    else
        k = ifftshift(-n/2:n/2-1);
    end
    
    % Spectral multipliers
    m = (k.^2)';
    
    if isempty(M);
        M = m;
    else
        % Create the permutation vector
        M_perm = circshift(1:length(size(M))+1,[0 1]);
        % Permute a singleton dimension to the beginning of M
        M = permute(M,M_perm);
        % Add m along the first dimension of M
        M = bsxfun(@plus,M,m);
    end
end
% Reverse M to match X (since new dimensions were being added form the left)
M = permute(M,length(size(M)):-1:1);
assert(isequal(size(M),size(X)),'Bad programming error'); % This should never occur
 
% Shape the amplitude multipliers by beta/4 which corresponds to shaping the power by beta
M = M.^(beta/4);
 
% Set the DC component to zero
M(1,1) = 0;
 
%% Multiply X by M in frequency domain
Xstd = std(X(:));
Xmean = mean(X(:));
X = real(ifftn(fftn(X).*M));
 
% Force zero mean unity standard deviation
X = X - mean(X(:));
X = X./std(X(:));
 
% Restore the standard deviation and mean from before the spectral shaping.
% This ensures the random sample from randn is truly random. After all, if
% the mean was always exactly zero it would not be all that random.
X = X + Xmean;
X = X.*Xstd;
end

Answer 1

Here is my solution:

My "min/max" code (lines 6 and 7) was bad. I wanted to divide all values in the matrix by the single largest value in the matrix so that all values would be between 0 and 1. Because I used max() improperly, I was stepping through the max value of each column and using that as my divisor; thus the vertical stripes.

In the end this is what my code looks like. X is the 3 dimensional matrix:

minVal = min(X,[],'all'); % Get the lowest value in the entire matrix
X = X - minVal; % Set min value to zero
maxVal = max(X,[],'all'); % Get the highest value in the entire matrix
X = X ./ maxVal; % Set max value to one