clg
% this gives a frequency that is exactly one of the FFT freqs
t1=[0:10*256-1]*2*pi/8;
a1=1;
x=a1*exp(i*t1);


sig=2;
x=x+sig*randn(size(x))+i*sig*randn(size(x));
subplot(211);
plot([1:256],real(x(1:256)));
v=axis;
axis([0 256 v(3) v(4)]);
title('Re(X) Time Series');
xlabel('index k');


axis;
subplot(212);
y=fft(x);
plot(2*log10(abs(y)));
v=axis;
axis([0 length(y) v(3) v(4)]);
title(['Periodogram of X with N=2560 and sigma_Y = ',num2str(sig*sqrt(2))]);
xlabel('Frequency index j');
ylabel('LOG_10(ABS(Y^2))');
nper=2560/8;
num = abs(y(nper+1));
y=[y(1:nper),y(nper+2:length(y))];
avnoise=mean(abs(y).*abs(y));
truepar = (length(x)*a1*a1)/(2*sig*sig);
% the 2 appears  in  truepar because there are indep real and imaj parts
estpar = (num^2)/avnoise;
disp(['truepar: ',num2str(truepar),' estpar: ',num2str(estpar)]);
%gtext(['LOG10 average of noise-only squares: ',num2str(avnoise),...
%' signal ',num2str(2*log10(abs(y(nper))))]);
%axis;


print -deps pgram5