空間周波数フィルタリング
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() の出力は複素数ですが、これの絶対値がマグニチュードで、マグニチュードの二乗がパワーです。