书上看到一段二维傅里叶变换的标准代码，总结如下。先来看看测试程序

function fft2_example
% fft2_example - 2D FFT example
clc
close all
clear all

% 产生矩孔
wx = 0.1;                    % rect x half-width (m)
wy = 0.05;                    % rect y half-width (m)
w = [wx wy];
L = 2;                        % side length x&y (m)
M = 200;                      % samples/side length

[g, x, y] = DiffractionPlane_Rect(L, M, w);

figure(1)
imagesc(x,y,g);             % image display
colormap('gray');           % linear gray display map
axis square;                % square figure
axis xy                     % y positive up
xlabel('x (m)'); ylabel('y (m)');

figure(2)
plot(x,g(M/2+1,:));         % x slice profile
xlabel('x (m)');
axis([-1,1,0,1.2]);

% 求二维傅里叶变换
G = FFT2D(g, L, L, 1);

% 求有效带宽(98%意义上的)
bandwidth = CalcEffBandWidth(G, 1/L, 1);
fprintf('有效带宽为:%gcyc/m\n', bandwidth);
end

%% 矩孔
function [u1, x1, y1] = DiffractionPlane_Rect(L1, M, w)
if length(w) == 1
    w(2) = w;
end

wx = w(1);                      % rect x half-width (m)
wy = w(2);                      % rect y half-width (m)

dx = L1/M;                      % src sample interval
x1 = -L1/2:dx:L1/2-dx;          % src coords
y1 = x1;
[X1, Y1] = meshgrid(x1, y1);

u1 = rect(X1/(2*wx)).*rect(Y1/(2*wy)); % src field
end

运行结果如下

前两幅图为空域图像及其中央剖面，中间三幅为其二维傅里叶变换，最后一幅意义如标题所示，从而可以求出原图像的有效带宽。

程序中具体使用的函数如下，已标准化

%% 矩孔
function [u1, x1, y1] = DiffractionPlane_Rect(L1, M, w)
if length(w) == 1
    w(2) = w;
end

wx = w(1);                      % rect x half-width (m)
wy = w(2);                      % rect y half-width (m)

dx = L1/M;                      % src sample interval
x1 = -L1/2:dx:L1/2-dx;          % src coords
y1 = x1;
[X1, Y1] = meshgrid(x1, y1);

u1 = rect(X1/(2*wx)).*rect(Y1/(2*wy)); % src field
end

function [G, fx, fy] = FFT2D(g, L, H, bshow)
% fcn: 求二维傅里叶变换
% Input: g - 二维图像；L，H - 二维图像的长宽；bshow - 控制是否显示频域图像
% Output: G - 变换结果；fx，fy - 频域的坐标
% Usage: 见fft2_example
% Date:  2011/8/26
% Author: GJ
% Email: guojiegj@mail.ustc.edu.cn
% Blog: http://blog.sina.com.cn/gjchunqiu

if nargin < 3
    error('GJ_Tips: Lack of parameters!');
elseif nargin == 3
    bshow = 0;
end

[M, N] = size(g);

dx = L/N;             % sample interval (m)
dy = H/M;

g0 = fftshift(g);     % shift
G0 = fft2(g0)*dx*dy;  % 2D fft and dxdy scaling
G = fftshift(G0);     % center

fx = -1/(2*dx):1/L:1/(2*dx)-(1/L); % x freq coords
fy = -1/(2*dy):1/H:1/(2*dy)-(1/H); % y freq coords

if bshow ~= 0
    figure
    surf(fx,fy,abs(G))  % display transform magnitude
    camlight left; lighting phong
    % colormap('gray')
    shading interp
    ylabel('fy (cyc/m)'); xlabel('fx (cyc/m)');
    title('magnitude');
   
    figure
    imagesc(fx, fy, nthroot(abs(G), 3));
    axis square; axis xy;
    colormap('gray');
    ylabel('fy (cyc/m)'); xlabel('fx (cyc/m)');
    title('magnitude');
   
    figure
    plot(fx,abs(G(M/2+1,:))); % fx slice profile
    title('magnitude');
    xlabel('fx (cyc/m)');
end
end

%% 求有效带宽(98%)
function [bandwidth, G2_power] = CalcEffBandWidth(G, dfx, bshow)
% fcn: 求二维傅里叶变换
% Input: G - 二维傅里叶变换结果；dfx - 频域单位（cyc/m）；bshow - 控制是否显示频域图像
% Output: bandwidth - 有效带宽; G2_power - 一定半径内的总能量
% Note: 要求G的横纵坐标fx，fy的采样间隔相同
% Usage: 见fft2_example
% Date:  2011/8/26
% Author: GJ
% Email: guojiegj@mail.ustc.edu.cn
% Blog: http://blog.sina.com.cn/gjchunqiu

G2 = abs(G).^2;             % spectral power
G2_sum = sum(G2(:));        % total power
G2_sum98 = G2_sum * 0.98;   % 98% effective bandwidth

[M, N] = size(G2);
Mhalf = floor((M+1)/2);
Nhalf = floor((N+1)/2);
y = -Mhalf:Mhalf-1;
x = -Nhalf:Nhalf-1;
[X, Y] = meshgrid(x, y);
R = sqrt(X.^2 + Y.^2);

rmax = floor(sqrt(Mhalf^2 + Nhalf^2)) + 1;
G2_power = zeros(rmax+1, 1);
for i = 0:rmax
    G2_i = G2.*(R<=i);
    G2_power(i+1) = sum(G2_i(:));
end

bandwidth = find(G2_power>=G2_sum98, 1, 'first') - 1;
bandwidth = bandwidth * dfx;

% show
if bshow ~= 0
    figure
    plot((0:rmax)*dfx, G2_power, '-o', ...
        [0 rmax]*dfx, [G2_sum98 G2_sum98], 'r', ...
        bandwidth, G2_sum98, '*r')
    xlabel('fr(cyc/m)');
    ylabel('power');
    title('随频域半径的增大，圆内包含能量的变化')
    text(0, G2_sum98*1.02, '98%')
end
end

function [out]=rect(x)

% fcn: rectangle function

out=abs(x)<=1/2;

end