clear all;
close all;

bits=1; % PN generator seed

Nt=1024; % Number of total data points 
Np=32; % Number of data points in each block of direct PSD estimate
     % 2Np-1 samples of total-data autocorrelation used in autocorr PSD estimate
     % Np must be element {2, 4, 8, ..., Nt/2}
Nb=round(Nt/Np); % Number of blocks
Na=2*Np-1; % Number of autocorrelation points kept
fftsize=Nt;

% Windows

zeropts=round(Nt-Np);
apsdwin=[zeros(1,zeropts) ones(1,Na) zeros(1,zeropts)];
dpsdwin=hamming(Np)';%rectwin(Np)';
fulldpsdwin=rectwin(round(0.5*Nt))';

% PN sequence generator

bitsarr=cell(1,Np*Nb); % PN generator bits in binary string form
B=16; % Number of bits in shift register
seq=zeros(1,Np*Nb);

% Generate "random" bits and sequence

for q=1:Nt
   nextbit=xor(bitand(bits,2^15),bitand(bits,2^1));
   nextbit=xor(bitand(bits,2^0),nextbit);
   bits=bitxor(bitshift(bits,1,B),nextbit);
   if nextbit==0
      seq(q)=1;
   else
      seq(q)=-1;
   end
   bitsarr{q}=dec2bin(bits,B);
end

disp('Some consecutive values of the shift register:');
disp(sprintf('%s\n',bitsarr{40:50}));

% Calculate psd of total PN sequence

pnpsd=fft(seq,fftsize); pnpsd=pnpsd.*conj(pnpsd);

% Calculate coloring filter and filter PN sequence

b=[32767 0 0 0 13421]/32767;
a=19345/32767;
h=freqz(a,b,fftsize,'whole');
seq=filter(a,b,seq);

% Estimates from Nt/2-point block of data

blk=seq(1:round(0.5*fftsize));

fulldpsd=fft(blk.*fulldpsdwin,fftsize);
fulldpsd=fulldpsd.*conj(fulldpsd);

fullapsd=abs(fft(xcorr(blk,blk,'biased'),fftsize));

% First-block-only estimates

blk=seq((Nt-Np+1):(Nt));

dpsd1=fft(blk.*dpsdwin,fftsize);
dpsd1=dpsd1.*conj(dpsd1);

apsd1=abs(fft(xcorr(blk,blk,'biased'),fftsize));

dpsdsum=zeros(1,fftsize);

% Estimates from total data (all blocks)

for q=1:Nb
   blk=seq([((q-1)*Np+1):q*Np]);
   dpsdsum=dpsdsum+abs(fft(blk.*dpsdwin,fftsize)).^2;
end
dpsd=dpsdsum/Nb;

apsd=abs(fft(xcorr(seq,seq,'biased').*apsdwin,fftsize));

% Plots

str=num2str(round(0.5*fftsize));
subplot(211); plot(fulldpsd*2/fftsize);
title(['Direct PSD, ' str '-point block']);
set(gca,'XLim',[0 Nt]);
subplot(212); plot(fullapsd);
title(['Autocorr PSD, ' str '-point block'])
set(gca,'XLim',[0 Nt]);

figure; str=num2str(Np);
subplot(211); plot(dpsd1/Np);
title(['Direct PSD, first ' str '-point block only']);
set(gca,'XLim',[0 Nt]);
subplot(212); plot(apsd1);
title(['Autocorr PSD, first ' str ...
       '-point block only, all autocorrelation points']);
set(gca,'XLim',[0 Nt]);

figure;
subplot(211); plot(dpsd/Np);
title(['Direct PSD, all ' str '-point blocks']);
set(gca,'XLim',[0 Nt]);
subplot(212); plot(apsd);
title(['Autocorr PSD, ' num2str(Na) ...
       ' autocorrelation points, entire sequence used']);
set(gca,'XLim',[0 Nt]);

figure; str=num2str(fftsize);
subplot(211); plot(pnpsd/fftsize);
title(['Spectrum of ' str '-point PN sequence']);
set(gca,'XLim',[0 Nt]);
subplot(212); plot(abs(h).^2);
title(['Impulse response of coloring filter']);
set(gca,'XLim',[0 Nt]);