clear;              %功率譜估計
Fs=10000;            %采樣頻率    
L=2;                %數據長度
n=0:1/Fs:L;         
N=length(n); 
x1=sin(2*pi*40*n);
e=randn(1,N);%白噪聲
x=x1+e;
k=100;
m=(N-1)/k;
xk=zeros(k,m);
X=zeros(k,2*m);
Rx=zeros(k,4*m-1);
M=zeros(k,2*m);
for i=1:k
    xk(i,:)=x((i-1)*m+1:i*m);
    z=zeros(1,m);
    X(i,:)=[xk(i,:) z];            %補零
    Rx(i,:)=xcorr(X(i,:),X(i,:));      %自相關 
    M(i,:)=abs(fft(Rx(i,:),2*m)); %做FFT并取模
end;
Sk=mean(M);
plot((1:2*m)*Fs/(2*m),Sk);  %把橫坐標轉換成頻率，并使數據一致