I wanted to take the convolution M different images (stored in an M x n x n) matrix with the same kernel in the Fourier domain, but I couldn’t figure out how to run fft2 on the second and third dimension of that matrix (computing a bunch of 2D FFTs in a single matlab call).

Turns out fft2 is a really simple function (type edit fft2 to see how it works) which involves calling fft on dimension 1 of an image, and then again on dimension 2. To take the 2d fft of all the images stored in the images matrix, then, proceed as follows:

imagesfft2 = fft(fft(images,[],2),[],3);

To convolve the kernel with the images, you need to multiply element wise with the 2d fft of the kernel. This is done cleanly using reshape and bsxfun:

kernfft2 = fft2(kern);
 kernfft2 = reshape(kernfft2,[1,size(kernfft2)])
 imagesfilt = bsxfun(@times,imagesfft2,kernfft2);

And to bring back to the image domain:

imagesconv = ifft(ifft(imagesfilt,[],2),[],3);

Easy breezy.