Fs=2000;
nFFT=1024;
n=0:1/Fs:1;
x=sin(2*pi*100*n)+4*sin(2*pi*500*n)+randn(size(n));
%直接法
X=fft(x,nFFT);
Pxx=abs(X).^2/length(n);
t=0:round(nFFT/2-1);
k=t/nFFT;
P=10*log(Pxx(t+1));
subplot(2,2,1);
plot(k,P);
%間接法
Cx=xcorr(x,'unbiased');
Cxk=fft(Cx,nFFT);
Pxx=abs(Cxk);
IP=10*log(Pxx(t+1));
subplot(2,2,2);
plot(k,IP);
%Bartlett法
K=4;
nFFTB=nFFT/K;
KK=1;
BPPx(nFFTB)=0;
while KK<=K,
    index=KK+(KK-1)*nFFTB;
    ii=1;
    while ii<=nFFTB,
        Bx(ii)=x(index+ii-1);
        ii=ii+1;
    end
    BX=fft(Bx,nFFTB);
    BPxx=abs(BX).^2/length(n);
    BPPx=BPPx+BPxx
    KK=KK+1;
end     
BPPx=BPPx/K;
bt=0:round(nFFTB/2-1);
bk=bt/nFFTB;
BP=10*log(BPPx(bt+1));
subplot(2,2,3);
plot(bk,BP);
%Welch法
M=64
L=(nFFT-M/2)/(M/2);
LK=1;
WPPx(M)=0;
while LK<=L,
    index=LK+(LK-1)*M/2;
    ii=1;
    while ii<=M,
        Wx(ii)=x(index+ii-1);
        ii=ii+1;
    end
    WX=fft(Wx,M);
    WPxx=abs(WX).^2/length(n);
    WPPx=WPPx+WPxx
    LK=LK+1;
end     
WPPx=WPPx/L;
wt=0:round(M/2-1);
wk=wt/M;
WP=10*log(WPPx(wt+1));
subplot(2,2,4);
plot(wk,WP);