randomly rough surface by ifft : real output from ifft

Question

I'm trying to generate a randomly rough isotropic surface with predefined roughness amplitudes (standard deviation of heights). Suppose I have the absolute values of fourier components of the surface in "H" and H has symmetry with respect to center of the array. Now, I need to define some random phase for the fourier series which are between -pi to pi (due to fftshift).

The question is how I can define these random numbers so that the final complex fourier components become complex conjugate? As in this case, the output of ifft will be real.

So far:

phi =  -pi + (pi+pi)*rand(n); % random phase 
[a,b] = pol2cart(phi,H);
  H_complex = complex(a,b); % the complex fourier transform composed of H & phi

  z = ifft2(ifftshift((H_complex)));

Answer 1

I found my answer, so maybe useful for someone else as well. I need to apply conjugate symmetry condition for 2D fourier transform, both in magnitude values and phase components. I applied it as below for a n*m matrix:

"H" is the magnitude:

H(1,1) = 0;
  H(1,m/2+1) = 0;
  H(n/2+1,m/2+1) = 0;
  H(n/2+1,1) = 0;
  H(2:end,2:m/2) = rot90(H(2:end,m/2+2:end),2);
  H(1,2:m/2) = rot90(H(1,m/2+2:end),2);
  H(n/2+2:end,1) = rot90(H(2:n/2,1),2);
  H(n/2+2:end,m/2+1) = rot90(H(2:n/2,m/2+1),2);

Similar operation but with negative sign applies to "phi" (phase):

phi(1,1) = 0;
  phi(1,m/2+1) = 0;
  phi(n/2+1,m/2+1) = 0;
  phi(n/2+1,1) = 0;
  phi(2:end,2:m/2) = -rot90(phi(2:end,m/2+2:end),2);
  phi(1,2:m/2) = -rot90(phi(1,m/2+2:end),2);
  phi(n/2+2:end,1) = -rot90(phi(2:n/2,1),2);
  phi(n/2+2:end,m/2+1) = -rot90(phi(2:n/2,m/2+1),2);

You can refer to this link for a good example of 2D fourier transform on a matrix and how the results should look like:

2D Fourier transform example