%直接法又稱為周期圖法，就是把x（n）的N個數據視為已能量有限的序列，直接計算
%x（n）的離散傅里葉變換，得X（k），然后再取其幅值的平方，作為功率譜的估計
t=linspace（011000）;
signal=4*sin（2*pi*50*t）+5*sin（2*pi*200*t）;
noise=randn（size（t））;
symbol=signal+noise;
Y=fft（symbol128）;
f=1000*（0:127）/128;
P1=Y.*conj（Y）/128;?%直接法
plot（f10*log10（P1（1:128）））;
xlabel（‘frequency‘）
ylabel（‘power‘）
title（‘直接法‘）
%或者用下面代碼
Fs=1000;%?采樣頻率
n=0:1/Fs:1;
xn=cos（2*pi*40*n）+3*cos（2*pi*100*n）+randn（size（n））;
window=boxcar（length（xn））;%矩形窗
nfft=1024;
[Pxxf]=periodogram（xnwindownfftFs）;
plot（f10*log10（Pxx））;

clear;
fs=500;??%采樣率
df=1;????%頻率分辨率
N=floor（fs/df）+1;%計算的序列點數
t=0:1/fs:（N-1）/fs;%截取信號的時間段
f=0:df:fs;%功率譜估計的頻率分辨率和范圍
xt=sin（2*pi*50*5）+2*sin（2*pi*130*t）+randn（1length（t））;
%截取時間段上的離散信號樣本序列



%利用周期圖法進行功率譜估計，但是其得出的功率譜很不光滑，相應的估計協方差比較大
%增加采樣點數也不能使周期圖變得更加平滑，這是周期圖的缺點，在后面對其改進。
Px=abs（fft（xt））.^2/（N^2）;%功率譜估計
Pav_tm=sum（xt.^2）/N;%在時域計算信號功率
Pav_fn=sum（Px）;%通過功率譜計算信號功率
figure（1）
subplot（221）???%作出功率譜密度圖
plot（f10*log10（Px））;
xlable（‘頻率（Hz）‘）;
Ylable（‘功率譜（dB）‘）;
title（‘周期圖得出的功率譜估計‘）;



%對周期圖進行改進的思想是將信號段進行估計，然后再將這些估計結果進行品均，
%從而減小估計的協方差，使估計功率譜圖變得平滑。
%本程序是將以上501點的信號分3段，分別作周期圖估計，然后平均。
Px=（abs（fft（xt（1:167）））.^2+abs（fft（xt（168:334）））.^2+abs（fft（xt（335:501）））.^2）/3/（（N/3）^2）;%分為三段
Pav_tm=sum（xt.^2）/N;%在時域計算信號功率
Pav_fn=sum（Px）;%通過功率譜計算信號功率
subplot（222）???%作出功率譜密度圖
plot（0:3:fs10*log10（Px））;
xlable（‘頻率（Hz）‘）;
Ylable（‘功率譜（dB）‘）;
title（‘分段估計得出的功率譜估計‘）;



%增加分段數可以進一步降低估計的協方差，然而每段中的數據點太少，就會使估計的
%頻率分辨率下降很多。在樣本信號序列總點數一定的條件下，可以采用使分段相互重疊
%的方法來增加分段數，從而保持每段中信號點數不變，這樣就在保證頻率分辨率的前提下
%進一步降低了估計