%Example done during Session 27 showing how autocorrelation
%sequence for sum of P sinusoids (SOS) can be extrapolated
%from the first P+1 values of the autocorrelation sequence
clf
clear
set(0,'defaultaxesfontsize',24);
%P=input('no. of complex sinewaves = ');
%omegas=input('frequencies of the sinewaves = ');
%amps=input('amplitudes of the sinewaves = ');
P=3;
omegas=[.45*pi .65*pi .75*pi];
amps=[.25 1 1];
rpos=zeros(512,1);
%
%create P+1 values of of the true autocorrelation sequence
%for R_{xx} (m) , m=0,1,2,...,P, based on SOS model
%
for i=1:P;
rpos(1:P+1,1)=rpos(1:P+1,1)+amps(1,i)^2*exp(j*omegas(1,i)*(0:P)');
end;
%
%solve for P linear prediction coefficients
%
R=toeplitz(conj(rpos(1:P,1)));
ar=-inv(R)*rpos(2:P+1,1);
%[U,S,V]=svd(R); Si=zeros(P,P);
%for im=1:P; Si(im,im)=1/S(im,im); end;
%ar=-U*Si*U'*rpos(2:P+1,1);
%
%extrapolate out to 512 values through linear prediction
%
for m=(P+2):512;
rpos(m,1)=0;
for k=1:P;
rpos(m,1)=rpos(m,1)-ar(k,1)*rpos(m-k,1);
end;
end;
%
% for different extrapolation lengths, compute DTFT of
%extrapolated autocorrelation sequence
%
Nlength=[P+1 32 64 128];
%plot extrapolated autocorrelation sequences
rpn=[conj(rpos((P+1):-1:2,1)) ; rpos(1:(P+1),1)];
s='kmcr';
plot((-P:P),abs(rpn),'k','Linewidth',4);
axis([-150 150 0 2.5])
title('Extrapolated Autocorrelation Sequence'),...
xlabel('lag value m'), ylabel('Magnitude'), grid
hold on
pause
for l=2:4;
r=[conj(rpos(Nlength(l):-1:2,1)) ; rpos(1:Nlength(l),1)];
x=(-Nlength(l)+1):(Nlength(l)-1);
plot(x,abs(r),s(l),'linewidth',4);
pause
end
hold off
Nlength=[P+1 32 64 128];
for l=1:4;
r=[conj(rpos(Nlength(l):-1:2,1)) ; rpos(1:Nlength(l),1)];
freqomega=linspace(-pi,pi,512);
fmag=20*log10(abs(fftshift(fft(r,512))));
fmag=fmag-max(fmag);
plot(freqomega,fmag,s(l),'LineWidth',4);
axis([-pi pi -30 0]);
title('Spectral Estimate'),...
xlabel('Omega'), ylabel('Magnitude'), grid
hold on
for p=1:P;
plot([omegas(1,p) omegas(1,p)],[-30 0],'b','LineWidth',3);
end
pause
end
clf
z=roots([1; ar]);
polar(0,1,'.')
hold on
plot(real(z),imag(z),'x','MarkerSize',18,'Linewidth',4);
plot(cos(omegas), sin(omegas),'o','MarkerSize',14,'Linewidth',4)
hold off


