clear

N = 512;

t = 0:(N-1);
fan1 = 0; fan2 = 0;
fs = 1000; %%%
f1 = 0.1*fs;
f2 = 0.4*fs;

x = sin(2*pi*f1/fs*t + fan1) + sin(2*pi*f2/fs*t + fan2);
y = awgn(x,6);

 plot(x); 
 figure,plot(y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% y 

y_z = fft(y,N);
%figure,plot(abs(y_z));
y_z = abs(y_z).^2/N;
%f1 = (0:fs/N:(N-1)*fs/N) - fs/2;
% figure,plot((y_z));
% %%figure,plot(f1,y);
% title('周期圖法');
%%axis([0 600 -4 130]);

figure,plot(10*log10(y_z));
title('周期圖法(dB)');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for m = 0:N-1
    x1 = y(1:N-m);
    x2 = y(1+m:N);
    r(m+1) = sum(x1.*x2)/N;
end;

a(1,1) = -r(2)/r(1);
ceta(1) = (1 - a(1,1).^2)*r(1);

k = 1;
while ceta(k) >= 0.5;
    k = k+1;
    temp = 0;
    for l = 1:(k-1)
        temp = temp + a(k-1,l)*r(k-l+1);
    end
    a(k,k) = -(r(k+1) + temp)/ceta(k-1);
    
    for i = 1:(k-1)
        a(k,i) = a(k-1,i) + a(k,k)*a(k-1,k-i);
    end
    
    ceta(k) = (1 - a(k,k)^2)*ceta(k-1);
end

ceta
k

N = N;
t = 0:(N-1);

w = 2*pi*t/N;

for n = 1:N
    temp = 0;
    for i = 1:k
        temp = temp + a(k,i)*exp(-w(n)*i*j);
    end
    yy(n) = 1/abs(1+temp);
end

% figure,plot((yy));
% title('levenson法');

figure,plot(10*log10(yy));
title('levenson法(dB)');