wf=pi/128;
x=ones(512);
wm=(0.5*(pi/6));
for n=1:512
    x(n)=cos(((pi/6)*n)+wm*(cos(wf*n)));
end
k=0.5;
wa=pi/128;
y1=ones(512);
for n=1:512
    y1(n)=(1+k.*cos(wa.*n)).*(x(n));
end
plot(y1)
for n=2:512
   y(n)=y1(n)-y1(n-1);
end
plot(y)
hold on
    for n=2:511
    ey(n) = (y(n)^2) - ((y(n-1)) * (y(n+1)));
end
plot(ey,'g')
hold on
for n=2:511
    ex(n) = (x(n)^2) - ((x(n-1)) * (x(n+1)));
end
plot(ex,'r')
hold on
for  n=2:510
   Q(n)=acos(1-(((ey(n)+ey(n+1))/(4*(ex(n))))));
end
plot(Q,'m')
hold on
for n=2:510
    j(n)=sqrt(ex(n)/(1-(1-(((ey(n)+ey(n+1))/(4*(ex(n))))))^2));
end
plot(j)