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

y = sign(rand(N, K)-0.5);
[f, p, pint] = avg_spectrum(y);

plot(f, 10*log10(p), f, 10*log10(pint));
xlabel('Frequency   [ f/f_s ]');
ylabel('Amplitude   [ dBWN ]');
axis([min(f) max(f) -20 10]);

title(sprintf('Random binary +/-1 sequence (%d averages)', K));