function imageFiltered = SpatialFilter(image, cutoff, pass)

% image: square array
% cutoff: cycles per image
% pass: 1) high-pass. 2) low-pass.

%% FFT
dc = mean(image(:));
image = image – dc;
imageFFT = fft2(image);

%% Power spectrum
powerS = log10(abs(fftshift(imageFFT).^2));
figure;
imagesc(powerS, [min(powerS(:)) max(powerS(:))]);
colormap(gray(256));
axis image;

%% Spatial Filter
[m, n] = size(image);
[x, y] = meshgrid(1:m, 1:n);
if pass == 1
filter = ((x-((m+1)/2)).^2 + (y-((n+1)/2)).^2) > cutoff^2;
else
filter = ((x-((m+1)/2)).^2 + (y-((n+1)/2)).^2) < (cutoff+1)^2;
end
filter = ifftshift(filter);

%% Filtering and inverse FFT
imageFilteredFFT = imageFFT.*filter;
imageFiltered = ifft2(imageFilteredFFT, ‘symmetric’);
imageFiltered = imageFiltered + dc;


関数は以上です。Psychtoolbox 関数は使用していません。フィルタをかける画像が必要なので、ここでは cameraman.tif を読み込んで、変数の型を double に変更します。
image = imread(‘cameraman.tif’);
image = cast(image, ‘double’);
imagesc(image, [0 255]);
colormap(gray);
axis image;

imageFiltered = SpatialFilter(image, 30, 1); を実行すると、配列変数 image に対して 30 cycles per image 未満の空間周波数をカットするハイパスフィルタをかけます（30 cycles はカットされません）。また、SpatialFilter() はパワースペクトルを表示します。

出力されたハイパス画像を表示するには以下を実行します。
imagesc(imageFiltered, [min(imageFiltered(:)) max(imageFiltered(:))]);
colormap(gray(256));
axis image;

元画像は256階調のグレースケールですが、フィルタをかけた画像の最小値は −31.7 、最大値は 301 になっていることに注意してください。実験刺激をフィルタリングする場合、フィルタされた画像の値が出力範囲（例えば0~255）におさまっている必要があります。

imageFiltered = SpatialFilter(image, 30, 2); を実行すると、配列変数 image に対して 30 cycles per image を超える空間周波数をカットするローパスフィルタをかけます（30 cycles はカットされません）。

2021年12月5日追記　マグニチュードスペクトルをプロットしていたので、パワースペクトルをプロットするように変更しました。fft() の出力は複素数ですが、これの絶対値がマグニチュードで、マグニチュードの二乗がパワーです。

function gaussian = GenerateGaussian(sd, imageSize)

% sd: pixel
% imageSize: pixel

[x, y] = meshgrid(1:imageSize, 1:imageSize);
x = x-(imageSize+1)/2; % cartesian coordinates
y = y-(imageSize+1)/2; % cartesian coordinates
radius = sqrt(x.*x + y.*y); % radial polar coordinates

gaussian = exp(-(radius.^2)/(2*sd^2)); % Gaussian: 0 to 1


関数は以上です。例えば、gaussian = GenerateGaussian(24, 200); を実行すると、200 x 200 の配列変数 gaussian に標準偏差 24 pixels のガウス関数の値が入ります。出力 gaussian の最小値は 0、最大値は 1 です。imshow(gaussian, [0 1]); を実行すると、作成したガウス関数を画像で確認できます。

Psychtoolbox 関数は使用していません。前回の GenerateGrating() と GenerateGaussian() を組み合わせるとガボアパターンを作成することができます。

function grating = GenerateGrating(wavelength, orientation, phase, imageSize)

% wavelength: pixel
% orientation: degree
% phase: degree
% imageSize: pixel

[x, y] = meshgrid(1:imageSize, 1:imageSize);
x = x-(imageSize+1)/2;
y = y-(imageSize+1)/2;

oriRad = -orientation/180*pi; % convert to radian
phaRad = phase*pi/180; % convert to radian
cpiRad = 2*pi/wavelength;

rotXY = (cos(oriRad).*x+sin(oriRad).*y); % rotate coordinates
grating = sin(rotXY*cpiRad+phaRad); % Sinusoidal grating: -1 to 1


関数は以上です。例えば、grating = GenerateGrating(40, 0, 90, 200); を実行すると、200 x 200 の配列変数 grating に波長 40 pixels で方位 0 度、位相 90 度のサイン関数の値が入ります。出力 grating の最小値は −1、最大値は 1 です。imshow(grating, [-1 1]); を実行すると、作成したグレーティングを画像で確認できます。

Psychtoolbox 関数は使用していませんが、imshow() を使うには image processing toolbox が必要です。imshow()を使えない場合は image(grating, ‘CDataMapping’, ‘scaled’); colormap(gray); axis image; で代用できます。

『はじめよう実験心理学（2015, 勁草書房）』の第３章でもグレーティングの作成について解説していますが、そちらは初心者向けのシンプルな関数となっています。