B = 1;
A = [1, -0.8];
tn = 0:1023;
wn = gaussnoise(0,1,1024); %generate the guauss distributed noise w(n)
xn = filter(B, A, wn);   %generate the signal x(n);

figure(1)
subplot(2,1,1)
plot(tn, wn)
title('噪聲信號w(n)時域波形圖');
ylabel('白噪聲w(n)')
xlabel('Time domain')
subplot(2,1,2)
plot(tn, xn)
title('信號x(n)時域波形圖');
ylabel('信號x(n)')
xlabel('Time domain')

%Pxx(exp(jw))=(abs(1/(1-0.8exp(-jw))))^2
Wn = linspace(0,pi,512);
Pxx=(1./(abs(1-0.8*exp(-j*Wn)))).^2;
figure(2)
subplot(3,1,1)
plot(Wn/pi, Pxx);
title('信號x(n)理想功率譜密度')
ylabel('信號功率譜Pxx(exp(jw))')
xlabel('歸一化頻率w(rad/s)')


subplot(3,1,2)
%figure(3)
[P0, Wn0] = periodogram_yuqd(xn, 1024);
plot(Wn0/pi,P0)
title('N＝1024時信號x(n)周期圖譜估計')
ylabel('信號功率譜Pxx(exp(jw))')
xlabel('歸一化頻率w(rad/s)')
%figure(4)
subplot(3,1,3)
[P1, Wn1] = periodogram_yuqd(xn(1:256), 1024);
plot(Wn1/pi,P1)
title('N＝256時信號x(n)周期圖譜估計')
ylabel('信號功率譜Pxx(exp(jw))')
xlabel('歸一化頻率w(rad/s)')

%Pxx(exp(jw))=(abs(1/(1-0.8exp(-jw))))^2
Wn = linspace(0,pi,512);
Pxx=(1./(abs(1-0.8*exp(-j*Wn)))).^2;
figure(3)
subplot(3,1,1)
plot(Wn/pi, Pxx);
title('信號x(n)理想功率譜密度')
ylabel('信號功率譜Pxx(exp(jw))')
xlabel('歸一化頻率w(rad/s)')
subplot(3,1,2)

[P0, Wn0] = periodogram_seg_yuqd(xn(1:256), 2, 1024);
plot(Wn0/pi,P0)
title('N＝256,L=2時信號x(n)平均周期圖譜估計')
ylabel('信號功率譜Pxx(exp(jw))')
xlabel('歸一化頻率w(rad/s)')

subplot(3,1,3)
[P1, Wn1] = periodogram_seg_yuqd(xn(1:256), 8, 1024);
plot(Wn1/pi,P1)
title('N＝256, L=8時信號x(n)平均周期圖譜估計')
ylabel('信號功率譜Pxx(exp(jw))')
xlabel('歸一化頻率w(rad/s)')