随机信号的功率谱密度估计和相关函数随机信号的功率谱密度估计和相关函数1．实验目的了解估计功率谱密度的几种方法，掌握功率谱密度估计在随机信号处理中的作用。
⒉实验原理随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。
功率谱密度简称为功率谱，是自相关函数的傅里叶变换。
对功率谱密度的估计又称功率谱估计。
1.线性估计法（有偏估计）：线性估计方法是有偏的谱估计方法，谱分辨率随数据长度的增加而提高。
包括自相关估计、自协方差法、周期图法。
2.非线性估计（无偏估计）：非线性估计方法大多是无偏的谱估计方法，可以获得高的谱分辨率。
包括最大似然法、最大熵法⒊实验任务与要求1. 所有功能均用matlab仿真。
2. 输入信号为：方波信号+n(t)，方波信号信号基频1KHz，幅值为1v，n(t)为白噪声。
3. 编写自相关估计法、自协方差法、周期图法、最大似然法、最大熵法的matlab 程序。
正确的运行程序。
4. 必须用图示法来表示仿真的结果。
对几种功率谱估计的方法进行比较分析，发现它们各自有什么特点？。
5. 按要求写实验报告。
4.Matlab程序如下：生成输入信号：clear;fs=1024;%设采样频率为1024n=0:1/fs:1;N=length(n);W=2000*pi;%因方波频率F=1000HZ所以角频率W=2000piX1n=square(W*n);%方波信号X2n=randn(1,N);%白噪声信号xn=X1n+X2n;%产生含有噪声的信号序列XNsubplot(3,1,1)plot(n,xn);xlabel('n')ylabel(‘输入信号’)%绘输入信号图（1）.周期图法：fs=4000;n=0:1/fs:1;N=length(n);W=2000*pi;x1n=square(W*n);x2n=randn(1,N);xn=x1n+x2n;subplot(3,1,1)plot(n,xn);Nfft=256;N=256;%傅里叶变换的采样点数256Pxx=abs(fft(xn,Nfft).^2)/N;f=(0:length(Pxx)-1)*fs/length(Pxx);subplot(3,1,2),plot(f,10*log10(Pxx)),%转成DB单位xlabel('频率/HZ'),ylabel('功率谱/db'),title('周期图法');（2）.相关函数法：fs=1000;n=0:1/fs:1;N=length(n);W=2000*pi;x1n=square(W*n);x2n=randn(1,N);xn=x1n+x2n;subplot(3,1,1)plot(n,xn);%输入信号m=-100:100[r,lag]=xcorr(xn,100,'biased')%求XN的自相关函数R，biased为有偏估计lag为R 的序列号subplot(3,1,2)hndl=stem(m,r);%绘制离散图，分布点从-100—+100set(hndl,'Marker','.')set(hndl,'MarkerSize',2);ylabel('自相关函数R(m)')%利用间接法计算功率谱k=0:1000;%取1000个点w=(pi/500)*k;M=k/500;X=r*(exp(-j*pi/500).^(m'*k));%对R求傅里叶变换magX=abs(X);subplot(3,1,3)plot(M,10*log10(magX));xlabel('功率谱的改进直接法估计')（3）.自协方差法：clear all;fs=1000;n=0:1/fs:3;P=2000*pi;y=square(P*n);xn=y+randn(size(n));%绘制信号波形subplot(211)plot(n,xn)xlabel('时间（s）')ylabel('幅度')title('y+randn(size(n))')ymax_xn=max(xn)+0.2;ymin_xn=min(xn)-0.2;axis([0 0.3 ymin_xn ymax_xn]) %使用协方差法估计序列功率谱p=floor(length(xn)/3)+1;nfft=1024;[xpsd,f]=pcov(xn,p,nfft,fs,'half'); %绘制功率谱估计pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001); subplot(2,1,2)plot(f,xpsd)title('基于协方差的功率谱估计') ylabel('功率谱估计（db）') xlabel('频率（HZ）')grid on;ymin=min(xpsd)-2;ymax=max(xpsd)+2;axis([0 fs/2 ymin ymax])（4）.最大熵法fs=4000;n=0:1/fs:1;N=length(n);W=2000*pi;x1n=square(W*n);x2n=randn(1,N);xn=x1n+x2n;subplot(3,1,1)plot(n,xn);Nfft=256;%分段长度256[Pxx,f]=pmem(xn,14,Nfft,fs);%调用最大熵函数pmem，滤波器阶数14 subplot(2,1,2),plot(f,10*log10(Pxx)),title(' 最大熵法，滤波器14'),xlabel('频率HZ'),ylabel('功率谱db');（5）.最大似然法：fs=1000;n=0:1/fs:1;N=length(n);W=2000*pi;x1n=square(W*n);x2n=randn(1,N);xn=x1n+x2n;subplot(3,1,1)plot(n,xn);%估计自相关函数m=-500:500;[r,lag]=xcorr(xn,500,'biased');R=[r(501) r(502) r(503) r(504);r(500) r(501) r(502) r(503);r(499) r(500) r(501) r(502);r(498) r(499) r(500) r(501)]; [V,D]=eig(R);V3=[V(1,3),V(2,3),V(3,3),V(4,3)].'; V3=[V(1,4),V(2,4),V(3,4),V(4,4)].'; p=0:3;wm=[0:0.002*pi:2*pi];B=[(exp(-j)).^(wm'*p)];A=B;%最小方差功率谱估计z=A*inv(R)*A';Z=diag(z');pmv=1./Z;subplot(2,1,2)plot(wm/pi,pmv);title('基于最大似然的功率谱估计') ylabel('功率谱幅度（db）') xlabel('角度频率w/pi')5.设计思想随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。
功率谱密度简称为功率谱，是自相关函数的傅里叶变换。
对功率谱密度的估计又称功率谱估计。
平稳随机信号x(t)的（自）功率谱Sxx(ω)定义为(1)式中rxx(τ)为平稳随机信号的自相关函数。
对于离散情况，功率谱表示为(2)式中T为离散随机信号的抽样间隔时间。
当利用随机信号的N个抽样值来计算其自相关估值时，即可得到功率谱估计为(3)可见，随机信号的功率谱与自相关函数互为傅里叶变换的关系，这两个函数分别从频率域和时间域来表征随机信号的基本特征。
按上式计算功率谱估值，其运算量往往很大，通常采用快速傅里叶变换算法，以减少运算次数。
线性估计方法是有偏的谱估计方法，谱分辨率随数据长度的增加而提高。
非线性估计方法大多是无偏的谱估计方法，可以获得高的谱分辨率。
周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度N。
由于序列x(n)的离散傅里叶变换X(k)具有周期性，因而这种功率谱也具有周期性，常称为周期图。
周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。
因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。
由于随机起伏大，使用周期图不能得到比较稳定的估值。
自相关函数是描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度；互相关函数给出了在频域内两个信号是否相关的一个判断指标，把两测点之间信号的互谱与各自的自谱联系了起来。
它能用来确定输出信号有多大程度来自输入信号，对修正测量中接入噪声源而产生的误差非常有效。
维纳-辛钦定理：随机信号的相关函数与其功率谱是一傅立叶变换对,即相关函数的傅立叶变换为功率谱,而功率谱的逆傅立叶变换为相关函数。
最大似然法原理是让信号通过一个滤波器，选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过，同时，使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。
在这个条件下，信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。
可以证明,如果信号x是由一个确定性信号S 加上一个高斯白噪声n所组成，则上述滤波器的输出是信号S的最大似然估值。
如果n不是高斯噪声，则上述滤波器的输出是信号S的最小方差的线性的无偏估值。
最大熵法主要思想是，在只掌握关于未知分布的部分知识时，应该选取符合这些知识但熵值最大的概率分布。
因为在这种情况下，符合已知知识的概率分布可能不止一个。
熵定义的实际上是一个随机变量的不确定性，熵最大的时候，说明随机变量最不确定，也就是随机变量最随机，对其行为做准确预测最困难。
从这个意义上讲，最大熵原理的实质就是，在已知部分知识的前提下，关于未知分布最合理的推断就是符合已知知识最不确定或最随机的推断，这是我们可以作出的唯一不偏不倚的选择，任何其它的选择都意味着我们增加了其它的约束和假设。
6.实验过程中遇到的问题及解决（1）由于采样间隔Ts=1/fs太小，故所产生信号波形难于观察，所以采用fs=1024;%设采样频率为1024。
这样产生的波形便易于观察。
（2）首次实验时，设置的fs太高，可能由于电脑配置的问题无法运行，降低了频率后，问题得以解决。
（3）对于matlab中的函数使用不熟练，不了解其中的含义及其功能。