N  = 2^10;
K  = 30;                % number of averages;

sp = zeros(N/2, 1);
K = 30;
for k=1:K
    Y = sign(rand(N, 1)-0.5);
    s = abs(fft(Y/sqrt(N)));
    sp = sp + s(1:N/2);
    sum(s.^2)/N
end
sp = sp / K;
    
f = (0:length(sp)-1) / N;

plot(f, 20*log10(sp));
xlabel('Frequency   [ f/f_s ]');
ylabel('Amplitude   [ dBWN ]');
axis([min(f) max(f) -80 20]);
fixfig;
