Here is a vectorized code that doesn't have to be modified if you change M or N:
th=pi/180;
u=[40*th 50*th 60*th 70*th];
b=u;
K=length(u);
% Note that N<M coprime No
N=3;
M=4;
NE = N+M-1; % Total
distance = [0; reshape((1:floor(NE/2))*pi.*[N; M],[],1)];
distance(NE+1:end) = [];
Sto = exp(-1i*distance.*cos(u));
Ste = exp(-1i*distance.*cos(b));
xx = abs(sum(Sto,2));
yy = abs(sum(Ste,2));
Note that you have b=u above, so that Sto == Ste, xx == yy, correlation is 1, and mse is 0. I suspect that your mse calculation is not what you really want (see comments in the code below, and maybe test with a case where b ~= u, so that mse is not zero, e.g., b = u+1, to see the difference between the two mse calculations).
%%%%%%%%%%%%%%%%%%%%%%%%
% MSE
%%%%%%%%%%%%%%%%%%%%%%%%
correlation = corr(xx,yy)
% mse1 = zeros(1,NE);
% for pp=1:NE
% mse1=(xx(pp)-yy(pp)); % I assume you meant "mse1(pp)" here instead of "mse1"
% end
mse = sum(xx-yy,1).^2/NE % this is equivalent to what you had (with the mse1(pp) correction in the loop above) ...
mse = mean((xx-yy).^2,1) % ... but this is MSE
mse_corr = mse + correlation-1 % note that correlation is a real scalar, so norm(correlation) == correlation (i.e., norm() does nothing)
