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