Considering A and B your two matrices, you could use the following code:
%reshape each matrix
SIZE = size(A);
A_2d = reshape(A,[SIZE(1)*SIZE(2),SIZE(3)]);
B_2d = reshape(B,[SIZE(1)*SIZE(2),SIZE(3)]);
% Calculate Pearson's correlation. It can be easily computed by transoforming data
% to zero mean and unit standard deviation (i.e., zscore). Then it is just a dot product.
A_2d_z = zscore(A_2d');
B_2d_z = zscore(B_2d');
correlation_1d = sum(A_2d_z.*B_2d_z)/(SIZE(3)-1);
%reshape back to the original size
correlation_2d = reshape(correlation_1d,[SIZE(1),SIZE(2)]);
