To calculate the standard deviation, the mean must be calculated again. Try to combine this:
x2 = g_p(d-w+q-1:d+q-2,:);
mean_x = sum(x2, 1) / w;
xc = x2 - mean_x; % Auto-expand: >= R2016b
% xc = bsxfun(@minus, x2, mean_x);
std_x = vecnorm(xc) / sqrt(s - 1); % vecnorm: >= R2017b
% std_x = sqrt(sum(xc .* xc, 1)) / sqrt(s - 1);
for the mean only the first and the last element changed between the iterations. Use this detail:
mean_x = sum(g_p(d-w:d-1, :) / w; % For q=1
for q = 1:n1-d
...
mean_x = mean_x - (g_p(d-w+q-1, :) + g_p(d+q-1, :)) / w;
end
