功率谱密度的估计
原始波=余弦波+白噪声
这个实验采用了两个输入，一个是白噪声，一个是有用信号和噪声信号作为输入时，他们的功率谱密度的仿真图像，并将他们进行对比。

平稳随机信号的功率谱密度（PSD ）是相关序列的离散傅里叶变换：
()()jwm XX x P w r m e ∞
--∞=∑
采用间接法计算噪声信号的功率谱。

间接法，又称自相关法或者BT 法，在1985年由布莱克曼与图基首先开拓。

间接法的理论基础是维纳-辛钦定理。

他是由N 个观察值x(0),x(1),……,x(N-1),估计出自相关函数R （m ）,然后再求R （m ）的傅里叶变换作为功率谱密度的估计。

()(),||1M jw
jwm N m M S e R m e M N -=-=<=-∑
clear all;
randn('state',0)
NFFT=1024; %采样点数
Fs=1000; %取样频率（单位为Hz ） t=0:1/Fs:.2;
y1=cos(t*20*pi); %余弦序列
figure(1)
plot(t,y1);
ylabel('余弦序列');
grid on;
%余弦序列的图像：
%白噪声
m=(0:NFFT-1)/Fs;
y=0.1*randn(size(m)); %产生高斯白噪声。

figure(2);
plot(m,y);
title('白噪声波形');
grid on;
%白噪声的自相关函数
[cory,lags]=xcorr(y,200,'unbiased'); %计算白噪声的自相关函数
figure(3)
plot(lags,cory); %自相关函数(无偏差的)，其中，cory为要求的自相关函数，lag为自相关函数的长度。

title('白噪声相关函数');
grid on;
%白噪声的频谱
f=fft(cory);
k=abs(f);
fl=(0:length(k)-1)*Fs/length(k); %f1为他的横坐标，单位为Hz.
figure(4)
plot(fl,k);
grid on;
title('白噪声功率谱'); % 自相关函数的傅里叶变换：即功率谱密度。

%输入信号为余弦信号和噪声信号时的功率谱密度
x= cos(t*2*pi*200)+0.1*randn(size(t)); %输入的一个噪声和有用信号（此时设输入信号为此，更易观察图像）。

figure(5);
plot(t,x);
ylabel('输入序列');
grid on;
title('输入序列');
[Cx,y]=xcorr(x,'unbiased'); %计算序列的自相关系数（无偏估计）
figure(6) %自相关函数
plot(y,Cx);
ylabel('输入序列的自相关函数');
title('输入信号的自相关函数');
grid on;
X=fft(Cx,NFFT); %计算傅里叶变换 pxx=abs(X); %求解功率谱密度 n=0:round(NFFT/2-1); %round即位近似，四舍五入。

k=n*Fs/NFFT; %横坐标
p=10*log(pxx(n+1)); %纵坐标为对数功率谱密度（更为精确）。

figure(7)
plot(k,p);
xlabel('频率(kHz)');
ylabel('相对功率谱密度(dB/Hz)');
grid on;
title('信号的功率谱');
△一：下图为NFFT=1024不同Fs对应的图像：
由上可知理想信号的情况下，信号的频谱是光滑的，但是当加入噪声信号以后，图像变的不光滑。

所以在实际应用中，应该尽量减小噪声的影响。

５结论
在经典谱估计中，无论是间接法法还是其改进的方法，都存在着频率分辨率低、方差性能不好的问题，原因是谱估计时需要对数据加窗截断，用有限个数据或其自相关函数来估计无限个数据的功率谱，这其实是假定了窗以外的数据或自相关函数全为零，这种假定是不符合实际的，正是由于这些不符合实际的假设造成了经典谱估计分辨率较差。

另外，经典谱估计的功率谱定义中既无求均值运算又无求极限运算，因而使得谱估计的方差性能较差，当数据很短时，这个问题更为突出。

如何选取最佳窗函数、提高频谱分辨率，如何在短数据情况下提高信号谱估计质量，还需要进一步研究。

参考文献：《MA TLAB 语言高级编程》机械工业出版社
《MA TLAB7.x程序设计语言》第二版
《线性代数》。