It's because
randi(1,Nb)
generates an Nb*Nb array. You need
randi(1,Nb,1)
as below
clear;clc;
fc = 3e8;
Nb = 1000; %% Number of snapshots
c=3e8;
wavelength = 3e8/fc;
d = 0.5*wavelength;
N = 10;
beta=2*pi/wavelength; % %Wavenumber
A = 1;% %Signal amplitude
snr = 5;% %SNR(dB)
theta = [30 80];
sigma = sqrt((A^2)/(2*10^(snr/10)));% %Variance of noise
M = length(theta);% %Number of signals
a=(d)/(2*sin(pi/N));% %radius of the circular array
%source signal
for k=1:M
D(k,:) = randi(1,Nb,1);
S(k,:) =A*(2*D(k,:) - 1);
end
for i=1:M
for k=1:N;
phi(k)=(2*pi*(k-1))/N;
SteeringVector(k,i)=exp(j*(beta*a*cos(theta(i)*pi/180)*cos(phi(k))-sin(theta(i)*pi/180)*cos(phi(k))));
end
end
% White Gaussien noise
B = (sigma^2)*(randn(N,Nb)+j*randn(N,Nb))/sqrt(2);
%Array output:signal plus noise
X = SteeringVector*S+B;
% Estimation of the spatial correlation matrix of the observed signal
Rxx = X*X'/Nb;
theta1=[0:180];
for i=1:length(theta1)
A1=zeros([N 1]);
for k=1:N
phiii(k)=(2*pi*(k-1))/N;
A1(k,1)= exp(j*(beta*a*cos(theta1(i)*pi/180)*cos(phiii(k))-sin(theta1(i)*pi/180)*cos(phi(k))));
end;
PBeamforming (i)= real(diag(A1'*Rxx*A1))/(N^2);
end;
figure(1);
plot(theta1,10*log10(PBeamforming ));
title('Beamforming spectrum');
xlabel('Angle [degree]');
ylabel('PBeamforming [dB]');
grid on;
