Make ifft2 to give real output

Mona Mahboob Kanafi (view profile)

on 20 May 2015
Latest activity Commented on by Ralf Mackenbach

Ralf Mackenbach (view profile)

on 10 Apr 2019
Accepted Answer by Mona Mahboob Kanafi

Mona Mahboob Kanafi (view profile)

I'm trying to generate a randomly rough surface(isotropic surface). 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) as in:
phi = -pi + (pi+pi)*rand(n); % random phase
Now, I don't know what to do to make the output of my surface real..
Here is the gist of the code I used so far:
[a,b] = pol2cart(phi,H);
H_complex = complex(a,b); % the complex fourier transform composed of H & phi
z = ifft2(ifftshift((H_complex)));
If I need to apply conjugate symmetry to the fourier components, how should I do it? The output is complex and I need a real surface.

Mona Mahboob Kanafi (view profile)

on 21 May 2015
Edited by Mona Mahboob Kanafi

Mona Mahboob Kanafi (view profile)

on 21 May 2015

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: http://fourier.eng.hmc.edu/e101/lectures/image_processing/node6.html

Ralf Mackenbach

Ralf Mackenbach (view profile)

on 10 Apr 2019
Small correction to this. The condition that
H(1,1) = 0;
H(1,m/2+1) = 0;
H(n/2+1,m/2+1) = 0;
H(n/2+1,1) = 0;
Does not need to hold necessarily. So a better (examplary) version would be
% Choose even m and n
m = 100 ;
n = 100 ;
% Choose average value of surface
z0 = 1 ;
% Fill random examplary matrices
H = rand(n,m) ;
phi = 2*pi*rand(n,m) ;
% Apply symmetry to matrices
H(1,1) = z0*m*n ;
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);
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);
% Convert to polar
[a,b] = pol2cart(phi,H);
H_complex = complex(a,b);
% Inverse and plot
S = ifft2(H_complex, 'symmetric') ;
surf(S)

Thomas Koelen (view profile)

on 20 May 2015

use abs() or real()? depending on if you want only the real number or the absolute value.