%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