This can easily be done in 3 lines of code.
Line #1
First use im2col to reshape each distinct block neighbourhood of 20 x 20 into a single column. As such, the output of this will be a 400 x N matrix, where each column denotes a unique block neighbourhood that has been reshaped into a column. Each column will have 400 rows, as each neighbourhood has 400 elements (20 x 20). N would be the total number of unique blocks we have in your 40 x 100 image. This would amount to 10, as we can fit 2 blocks horizontally and 5 blocks vertically given the 20 x 20 requirement.
Line #2
What is great about the output of im2col is that the ith row of im2col tells you the ith element for every block in your image. As such, all you have to do next is take each row and average over all of the columns. The output will be a 400 x 1 vector that denotes your average FFT for all of the blocks. This can be achieved using mean and specifying that we want to average over the second dimension (second parameter is 2), which is the columns.
Line #3
We then need to reshape this back into a 20 x 20 matrix, so use reshape to do this. We specify that the output matrix is 20 x 20, given the 400 x 1 element vector.
One question that you may ask is whether or not this re-ordering is guaranteed to reorder our FFT block correctly. This is guaranteed because when im2col constructed each block into a column, it progresses in a column-major order. This means that for one column of blocks, we construct them on a row-by-row basis. Once we get our 20 x 20 set of distinct blocks, these blocks are arranged so that the are sampled in column major order. This means that a single 20 x 20 block gets constructed into a 400 x 1 column vector, where columns of the 20 x 20 block are stacked on top of each other from left to right. Therefore, by doing mean and reshape, the spatial locations for each block do correspond to each other and will thus produce the right answer.
Without further ado, here's the code:
colBlocks = im2col(ffts, [20 20], 'distinct'); %// Line 1
meanCol = mean(colBlocks, 2); %// Line 2
fftBlockAverage = reshape(meanCol, [20 20]); %// Line 3
Minor Side Effect
Because the FFT is complex valued in nature, by doing the average, you would average the real and imaginary components separately. This is how MATLAB handles the average of complex valued data. I'm not sure what analysis you'll be performing after you calculate the average 2D FFT block, but bear this in mind before you proceed any further with your analysis.
Sidenote
Divakar in an earlier answer created a more efficient implementation of im2col. This is especially useful if you don't have the Image Processing Toolbox installed. You can check out that implementation here. It has been shown that the timing between this function and MATLAB's im2col are magnitudes faster.
Benchmarking
As a bonus, here is a benchmark using his code. Timing results were taken using a 40 x 100 matrix where the im2col built-in function was timed, and Divakar's custom function after. Results show that his method is faster. This may be very useful when considering larger size images. However, if you are looking for succinctness, use what I have written. If you want something fast, use his method.
Benchmarking Code
%// Input Parameters
nrows = 20;
ncols = 20;
A = rand(40,100);
disp('------------------------- With IM2COL');
tic
B1 = im2col(A,[nrows ncols],'distinct');
toc,clear B1
disp('----------------- With CUSTOM-BUILT IM2COL');
tic
B2 = im2col_distinct(A,[nrows ncols]);
toc,clear B2
Results
------------------------- With IM2COL
Elapsed time is 0.026914 seconds.
----------------- With CUSTOM-BUILT IM2COL
Elapsed time is 0.004186 seconds.
