% example_fft2.m

%

% This MATLAB script file demonstrates the

% computation of 2D Fourier transforms.

%

% K. Beers. MIT ChE. 8/19/2002

% First, declare the grid.

% x-direction

Px = pi;

Nx = 2^6;

dx = 2*Px/Nx;

x_val = [0 : dx : (2*Px - dx)];

% y-direction

Py = pi;

Ny = 2^6;

dy = 2*Py/Ny;

y_val = [0 : dy : (2*Py - dy)];

% use meshgrid to produce X and Y

% values at 2-D grid points.  Here

% the (i,j) element of the matrix

% captures the value of a function

% at the position (x(i),y(j))

[X,Y] = meshgrid(x_val,y_val);

% We now create a pattern with a simple

% Fourier spectrum.

F = zeros(size(X));

F = F + cos(X) + cos(3.*X) + ...

    2*sin(X).*cos(2.*Y) + ...

    cos(4.*Y);

fun_name = 'cos(x) + cos(3x) + 2*sin(x)*cos(2y) + cos(4y)';

% make a plot of the function

figure;

subplot(2,1,1);

surf(X,Y,F);

xlabel('x'); ylabel('y'); zlabel('f(x,y');

title(['f(x,y) = ', fun_name]);

subplot(2,1,2);

contourf(X,Y,F,15);

colorbar;

xlabel('x'); ylabel('y'); zlabel('f(x,y');

% Now, generate the 2-D grid of sampled q.

% x-direction

dqx = pi/Px;

qx_crit = pi/dx;

qx_val = [0 : dqx : (2*qx_crit-dqx)];

% y-direction

dqy = pi/Py;

qy_crit = pi/dy;

qy_val = [0 : dqy : (2*qy_crit-dqy)];

[QX,QY] = meshgrid(qx_val,qy_val);

% We now invoke the discrete 2D FFT routine.

F_FT = fft2(F);

F_PS = real(conj(F_FT).*F_FT);

% Plot the Fourier transform

figure;

surf(QX,QY,F_PS);

xlabel('q_x'); ylabel('q_y');

zlabel('|F(q)|^2');

title(['Power spectrum of f(x,y) = ', ...

        fun_name]);

% We now make contour plots of the original

% function, the Fourier transform, and

% the mode of largest magnitude in the

% power spectrum.

figure;

contourf(QX,QY,F_PS,15);

axis([qx_val(1) 5 qy_val(1) 5]);

colorbar;

xlabel('q_x'); ylabel('q_y');

zlabel('|F(q)|^2');

title(['Power spectrum of f(x,y) = ', ...

        fun_name]);