随机信号的功率谱密度
随机信号的功率谱密度
估计和相关函数
随机信号的功率谱密度估计和相关函数
1．实验目的
了解估计功率谱密度的几种方法，掌握功率谱密度估计在随机信号处理中的作用。
⒉ 实验原理
随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。功率谱密度简称为功率谱，是自相关函数的傅里叶变换。对功率谱密度的估计又称功率谱估计。
1.线性估计法（有偏估计）：线性估计方法是有偏的谱估计方法，谱分辨率随数据长度的增加而提高。包括自相关估计、自协方差法、周期图法。 2.非线性估计（无偏估计）：非线性估计方法大多是无偏的谱估计方法，可以获得高的谱分辨率。包括最大似然法、最大熵法
⒊ 实验任务与要求
1. 所有功能均用matlab仿真。
2. 输入信号为：方波信号+n(t)，方波信号信号基频1KHz，幅值为1v，n(t)为白噪声。
3. 编写自相关估计法、自协方差法、周期图法、最大似然法、最大熵法的matlab程序。正确的运行程序。
4. 必须用图示法来表示仿真的结果。对几种功率谱估计的方法进行比较分析，发现它们各自有什么特点？。 5. 按要求写实验报告。
4.Matlab程序如下：
生成输入信号： clear;
fs=1024;%设采样频率为1024 n=0:1/fs:1; N=length(n);
W=2000*pi;%因方波频率F=1000HZ所以角频率W=2000pi X1n=square(W*n);%方波信号 X2n=randn(1,N);%白噪声信号 xn=X1n+X2n;
%产生含有噪声的信号序列XN
subplot(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;%傅里叶变换的采样点数256 Pxx=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—+100 set(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中的函数使用不熟练，不了解其中的含义及其功能。起初将floor(x)函数的功能与四舍五入混淆在一起，导致程序无法运行。但经过查阅资料了解到其功能是“下取整”，或者说“向下舍入”，即取不大于x的最大整数。 （4）起初不知道调用最大熵函数是pmem，使得程序无法继续编写，找到函数库后，解决了问题。
（5）自相关函数是描述随机信号X(t)在任意两个不同时刻t1，t2的取值之间的相关程度.设原函数是f(t)，则自相关函数定义为R(u)=f(t)*f(-t)，其中*表示卷积。自相关函数是信号间隔的函数，间隔有正负间隔，所以n个长度的信号，有2n-1
个自相关函数值，分别描述的是不同信号间隔的相似程度。
使用规则：C=XCORR(A,‟flag‟) 该函数返回长度为(2N-1)的自相关序列； 如果“flag”为“none”时，则计算序列的非归一化行相关；
如果“flag”为“biased”时，则计算自相关函数有偏估计；
如果“flag”为“unbiased”时，则计算自相关函数无偏估计；
如果“flag”为“coeff”时，则计算序列归一化行处理，使得对零滞后的样本其自相关序列恒为1。
7．结论
比较：周期图法与序列的频谱有对应的关系，能采用FFT算法快速实现。但在直接法功率谱估计中，要对无限长的序列加以矩形窗，这也意味着对自相关函数加窗，因此会产生频谱泄漏，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱模糊和失真，使周期图功率谱的分辨率较低。而最大熵法最大限度第保留了数据阶段后‘窗口’以外的信息，使估计谱的熵最大。最大似然估计的分辨率相对比较高。自相关法中X(n)的功率谱估计为P(k)=∑R(m)W.^(m*K)。当M较小时，计算量不算很大。
8. 心得与建议
心得：随机信号的功率谱密度估计和相关函数这个实验主要是了解估计功率谱密度的几种方法，掌握功率谱密度估计在随机信号处理中的作用。随机信号的功率谱密度用来描述信号的能量特征随频率的变化关系。功率谱密度简称为功率谱，是自相关函数的傅里叶变换。对功率谱密度的估计又称功率谱估计。最初看到随机实验选题的时候，对于实验的目的，实验的步骤在脑中没有概念。对于线性估计法和非线性估计法的差异没有明确的区分。通过查阅书籍，了解了五种估计方法的基本思想。但对于用matlab进行仿真估计，也十分不熟悉。于是，我们小组成员找出课本复习了matlab中的一些重点。通过我们不屑的努力，不断的修改程序，调试程序，终于得出了正确的仿真图形。通过这次试验，对功率谱密度有了更深刻更全面的认识，并且对matlab一些功能的使用也熟悉了许多。实践了理论，对理论有了更具体的认识。但仍看到了一些不足，理论基础不够扎实，对于matlab运用不够熟练。并且希望可以更深入的研究如何可以更加精确的描述功率谱密度的方法。
建议：通过完成随机信号分析实验这样在时间和内容上都很自由的开放性实验，锻炼了我们自己设计实验，完成实验的动手能力，同样也培养了我们合作的团队精神。但由于形式上的开放，也使我们在实验中遇到了许许多多的困难，起初在面对困难时觉得无从下手，但经过成员的商量，讨论，最终达成共识。当我们克服重重困难之后，发现这样形式的实验，使我们对于理论有了更加深刻的认识，从而能够将理论熟练的运用到实践中去。希望这样的实验能够多方面的开展，以便于我们进一步提高动手能力以及通过实验提高理论联系实际的能力，从而更深刻的了解理论的概念。