clc
clear?all;
M=10;?????%數(shù)據(jù)組數(shù)
N=1024;???%數(shù)據(jù)點(diǎn)數(shù)
nfft=2*N;?%FFT變換為復(fù)序列
fs=1;?????%歸一化采樣頻率



%%產(chǎn)生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;
%%
%真實(shí)功率譜syms?ff;
a=[zeros（1N）10.981.920.940.92];?%模型系數(shù)
b=xcorr（alength（a）-1）;??%自相關(guān)
c=（abs（fft（b）））.^2;
Px=-10*log10（c）;
figure;
plot（Px（1:N））;????%真實(shí)功率譜
title（‘真實(shí)功率譜‘）;xlabel（‘f（未歸一化）‘）;ylabel（‘P（dB）‘）

%%
%周期圖法功率譜估計(jì)
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法功率譜估計(jì)
figure;
L=[256?128?64];???%數(shù)據(jù)段長(zhǎng)度取值
for?i=1:1:length（L）
sum_Pxx=0;
for?j=1:1:N/L（i）
????window=hamming（L（i））;%漢明窗
????noverlap=length（wi