Code covered by the BSD License  

Highlights from
Huygens-Fresnel integral approximation, free-form apertures and rough surfaces

image thumbnail

Huygens-Fresnel integral approximation, free-form apertures and rough surfaces

by

 

Simulate wave propagation through free-form apertures, or off rough surfaces.

generateSurfaceRough(obj)
function [surface, nGrid] = generateSurfaceRough(obj)
% generate a random, band-limited surface with a specified roughess
% create a FFT with magnitude 1 and random phase, filter in the frequency
% domain, inverse transform to get an arbitrary surface, scale linearly to
% reach the desired R_a

% map to convenient variables
m = obj.nGrid;

% satisfy constraints: zero frequency real, all other values point
% symmetric and conjugate complex: fZ(x+x0,y+yo)==conj(fZ(-x+x0,-y+y0));
g = @(x) floor((x+mod(x,2))/2)*2-1; % get next smaller odd number
m0 = g(m(1));
n0 = g(m(2));

if m0 < 3 || n0 < 3
    displayMessage(-1, [obj.name ' ERROR: surface dimensions are to small']);
end

% general layout of the fft2:
%
% rZ1  Z2  Z2 ... ... ... ... ... ~Z2 ~Z2
%  Z3  Z4 ... ...  Z4  Z5  Z6 ... ...  Z6
%  .   .   .   .   .   .   .   .   .   .
%  .   Z4 ... ...  Z4  Z5  Z6 ... ...  Z6
%  .   Z7 ... ...  Z7 rZ8 ~Z7 ... ... ~Z7
%  .  ~Z6 ... ... ~Z6 ~Z5 ~Z4 ... ... ~Z4
%  .   .   .   .   .   .   .   .   .   .
% ~Z3 ~Z6 ... ... ~Z6 ~Z5 ~Z4 ... ... ~Z4
%
% where 'r' means real, and '~' means conjugate complex
% Z4 and Z6 are point symmetric around rZ8 

Z1 = 1;
Z2 = exp(1i * random('uniform', 0, 2*pi, 1, floor(n0/2)));
Z3 = exp(1i * random('uniform', 0, 2*pi, 1, floor(m0/2)))';
Z4 = exp(1i * random('uniform', 0, 2*pi, floor(m0/2), floor(n0/2)));
Z5 = exp(1i * random('uniform', 0, 2*pi, 1, floor(m0/2)))';
Z6 = exp(1i * random('uniform', 0, 2*pi, floor(m0/2), floor(n0/2)));
Z7 = exp(1i * random('uniform', 0, 2*pi, 1, floor(n0/2)));

fZ = [               Z4               Z5                Z6;...
                     Z7                1  conj(fliplr(Z7));...
      conj(rot90(Z6,2)) conj(flipud(Z5)) conj(rot90(Z4,2))];

% handle even numbered dimensions
if m(1) == m0       % m odd
    if m(2) ~= n0   % n even
        fZ = [ [Z3; 1; conj(flipud(Z3))], fZ];          % add first column     
    end
else                % m even
    fZ = [ Z2, 1, conj(fliplr(Z2)); fZ];                % add first row
    if m(2) ~= n0   % n even
        fZ = [ [Z1; Z3; 1; conj(flipud(Z3))], fZ];      % add first column 
    end
end

% apply frequency filter
[X Y] = ndgrid(-m(1)/2:(m(1)-1)/2, -m(2)/2:(m(2)-1)/2);
param(1,:) = obj.scaledFilterParameters(1,:) * m(1);
param(2,:) = obj.scaledFilterParameters(2,:) * m(2);
tp = obj.filter(X,Y,obj.filterParameters, param);

% transform back the artificial FFT2
global useGpu;
if useGpu && obj.useGpu
	if GPUisDoublePrecision
		displayMessage(2, [obj.name ': double precision IFFT (GPU)']);
		fZ_gpu = GPUdouble(ifftshift(fZ.*tp)); % try compute capabiltity 2.0
	else
		displayMessage(2, [obj.name ': single precision IFFT (GPU)']);
		fZ_gpu = GPUsingle(ifftshift(fZ.*tp)); % live with compute capability 1.x
	end
	Z_gpu = ifft2(fZ_gpu);
	Z = real(double(Z_gpu));
else
	displayMessage(2, [obj.name ': double precision IFFT (CPU)']);
	Z = real(ifft2(ifftshift(fZ.*tp)));
end

% subtract mean, and scale to reach 'ra'
Z = Z - mean(Z(:));
Zra = 1 / (m(1)*m(2)) * sum(abs(Z(:)));		% actual Ra
Zn = Z / Zra * obj.ra;

% get matrix size
nGrid = size(Zn);

% reshape Z
surface = zeros(3,numel(Zn));
surface(1,:) = reshape(X, 1, []) .* obj.resolution(1);
surface(2,:) = reshape(Y, 1, []) .* obj.resolution(2);
surface(3,:) = reshape(Zn, 1, []);

Contact us