clc
clear?all;
M=10;?????%數據組數
N=1024;???%數據點數
nfft=2*N;?%FFT變換為復序列
fs=1;?????%歸一化采樣頻率



%%產生AR（4）采樣樣本10*1024
for?n=1:M
??wn=wgn（1N-30）;??%白噪聲
for?i=1:4
??x（ni）=wn（i）;
end;
for?j=5:N
??x（nj）=-0.98*x（nj-1）-1.92*x（nj-2）-0.94*x（nj-3）-0.92*x（nj-4）+wn（j）;?
end;

end;
%%
%真實功率譜syms?ff;
a=[zeros（1N）10.981.920.940.92];?%模型系數
b=xcorr（alength（a）-1）;??%自相關
c=（abs（fft（b）））.^2;
Px=-10*log10（c）;
figure;
plot（Px（1:N））;????%真實功率譜
title（‘真實功率譜‘）;xlabel（‘f（未歸一化）‘）;ylabel（‘P（dB）‘）

%%
%周期圖法功率譜估計
figure;
sum_Pxx=0;
for?j=1:1:M
window=hamming（N）;%漢明窗
[Pxxf]=periodogram（x（j:）window2*Nfs）;?
sum_Pxx=sum_Pxx+abs（Pxx）;
end
f=f/max（f）/2;%歸一化
plot（f10*log10（sum_Pxx/M））;title（‘周期圖法‘）;xlabel（‘f（Hz）‘）;ylabel（‘P（dB）‘）;


%%
%WELCH法功率譜估計
figure;
L=[256?128?64];???%數據段長度取值
for?i=1:1:length（L）
sum_Pxx=0;
for?j=1:1:N/L（i）
????window=hamming（L（i））;%漢明窗
????noverlap=length（wi