%信號xn
N=1024;
n=[0:N-1];
f1=0.05;
f2=0.120;
wn=randn(1,N);
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn;

%周期圖法求功率譜
Pxx=10*log10(abs(fft(xn).^2)/N);
f=(0:length(Pxx)-1)/length(Pxx);
figure(1);
plot(f,Pxx);
xlabel('frequency');
ylabel('power spectrum(db)');
title('period');
grid;


%間接法求功率譜
r = zeros(N,1);        %求自相關(guān)函數(shù)
for m=0:16            %取m<<N-1
    for n=0:(N-1-m)
      r(m+1)=r(m+1)+xn(n+1)*xn(n+1+m);
    end
     r(m+1)=r(m+1)/N;
end
p=10*log10(fft(r,N));  %對自相關(guān)函數(shù)求傅立葉變換得功率譜
f=(0:length(p)-1)/length(p);
figure(2);
plot(f,p);
xlabel('frequency');
ylabel('power spectrum(db)');
title('間接法功率譜');
grid;

for mm=0:240            %取m=N-1
    for n=0:(N-1-mm)
      r(mm+1)=r(mm+1)+xn(n+1)*xn(n+1+mm);
    end
     r(mm+1)=r(mm+1)/N;
end
pm=10*log10(fft(r,N));  %對自相關(guān)函數(shù)求傅立葉變換得功率譜
f=(0:length(pm)-1)/length(pm);
figure(3);
plot(f,pm);
xlabel('frequency');
ylabel('power spectrum(db)');
title('間接法功率譜');
grid;


%bartlett法求功率譜
Ns=256;
w1=ones(1,Ns);%加矩形窗
pxx1=abs(fft(w1.*xn(1:256),Ns).^2)/norm(w1)^2;
pxx2=abs(fft(w1.*xn(257:512),Ns).^2)/norm(w1)^2;
pxx3=abs(fft(w1.*xn(513:768),Ns).^2)/norm(w1)^2;
pxx4=abs(fft(w1.*xn(769:1024),Ns).^2)/norm(w1)^2;
pxxt=10*log10((pxx1+pxx2+pxx3+pxx4)/4);
f=(0:length(pxxt)-1)/length(pxxt);
figure(4);
plot(f,pxxt);
xlabel('frequency');
ylabel('power spectrum(db)');
title('bartlett');
grid;


%welch法求功率譜
Ns=256;
w2=hanning(256)';%加漢寧窗
Pxx1=abs(fft(w2.*xn(1:256),Ns).^2)/Ns;
Pxx2=abs(fft(w2.*xn(129:384),Ns).^2)/norm(w2)^2;
Pxx3=abs(fft(w2.*xn(257:512),Ns).^2)/norm(w2)^2;
Pxx4=abs(fft(w2.*xn(385:640),Ns).^2)/norm(w2)^2;
Pxx5=abs(fft(w2.*xn(513:768),Ns).^2)/norm(w2)^2;
Pxx6=abs(fft(w2.*xn(641:896),Ns).^2)/norm(w2)^2;
Pxx7=abs(fft(w2.*xn(769:1024),Ns).^2)/norm(w2)^2;
Pxxt=10*log10((Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7)/7);
f=(0:length(Pxxt)-1)/length(Pxxt);
figure(5);
plot(f,Pxxt);
xlabel('frequency');
ylabel('power spectrum(db)');
title('welch');
grid;