%program p511
var=sqrt(1/exp(1.0));
v=var*randn(1,1000);

n=1:1000;
w0=1.4*pi;
xn=exp(j*pi*n-j*pi)+exp(j*w0*n-j*0.7*pi)+v;
m=-500:500;
[r,lag]=xcorr(xn,500,'biased');%corr

R=[r(501),r(502),r(503),r(504);
   r(500),r(501),r(502),r(503);
   r(499),r(500),r(501),r(502);
   r(498),r(499),r(500),r(501)];

[V,D]=eig(R);
V3=[V(1,3),V(2,3),V(3,3),V(4,3)].';
V4=[V(1,4),V(2,4),V(3,4),V(4,4)].';

p=0:3;
wm=[0:0.002*pi:2*pi];
B=[(exp(-j)).^(wm'*p)];
A=B.';

%p-mvdr
z=A'*inv(R)*A;
Z=diag(z');
pmv=1./Z;
plot(wm/pi,pmv);
xlabel('Angular frequency,ω/pi');
ylabel('Power Spectrum Magnitude (dB)');
title(' MVDR Spectral Estimate')