% matlab实现经典功率谱估计
%周期图法 直接法：
Fs=1000; %采样频率
n=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));%randn(size(n))返回和n同样维数的随机数组 
window=boxcar(length(xn)); %矩形窗 boxcar产生矩形框
nfft=1024;
[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 periodogram(xn,window,nfft,Fs)periodogram是用来计算功率谱密度的，参数中，X：所求功率谱密度的信号window：所使用的窗口，默认是boxcar，其长度必须与x的长度一致；nfft：采样点数；fs：采样频率。
plot(f,10*log10(Pxx));
% matlab实现经典功率谱估计
% 间接法：间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
clear all;
Fs=1000; %采样频率
n=0:1/Fs:1;
%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); 
nfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数
CXk=fft(cxn,nfft);
Pxx=abs(CXk);
index=0:round(nfft/2-1);
k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
plot(k,plot_Pxx);