Help increasing the speed of a piece of code

Dear all,   Plain summary: I have a temperature matrix and I need to find the average temperature of the 8 pixels surrounding each pixel.   I need to subdivide my 134x134 matrices into 3x3 sub matrices. For each submatrix I need to compute the mean and the standard deviation, excluding the central pixel from the calculation.   Then I need to find the values exceding the upper and lower boundaries, namely the mean + standard deviation and the mean - standard deviation, respectively. Then, I need to recalculate the mean excluding the values falling outside the lower and upper boundaries and assign the mean value to the central pixel location (in a new matrix).   The output will be a 134x134 matrix with every value (pixel) corresponding to the average of its surrounding pixels with the "outliers" excluded from the average.   A=magic(134); % Assuming this is my 134x134 Temperature matrix T_background = blockproc(A,[1,1],@BlkFun, 'BorderSize',[1,1],... 'TrimBorder',false, 'PadMethod',NaN); %where BlkFun= function v = BlkFun(s) px = s.data; if isequal(size(px),[3,3]) px(5) = []; av = nanmean(px); sd = nanstd(px); ix = px>=(av-sd) & px<=(av+sd); v = nanmean(px(ix)); else v = NaN; end end