随机过程基础(7)

互功率谱

类似从自相关函数推出互相关函数，也可以从单个随机过程的功率谱推出互功率谱。两个随机过程X(t)和Y(t)，它们的互功率谱定义是
$$G_{XY}(\omega )=E[\underset{T\to \infty }{\lim }\frac{1}{2T}{X}_T(\omega )Y_T^{\ast }(\omega )]$$
$$X_T(\omega )={\int }_{-T}^{T}X(t)e^{-j\omega t}dt$$
$$Y_T(\omega )={\int }_{-T}^{T}Y(t)e^{-j\omega t}dt$$
如果X(t)和Y(t)是联合平稳的，则$R_{XY}(\tau )$绝对可积。则它们的互功率谱和互相关函数是一对傅里叶变换对：
$$G_{XY}(\omega )={\int }_{-\infty }^{infty}{R}_{XY}(\tau )e^{-j\omega \tau }d\tau$$
$$R_{XY}(\tau )=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{G}_{XY}(\omega )e^{j\omega \tau }d\omega$$
互功率谱的性质：
$$G_{XY}(\omega )=G_{YX}^{\ast }(\omega )$$
$$|G_{XY}(\omega ){|}^2\le G_X(\omega )G_Y(\omega )$$

非平稳过程的功率谱

广义功率谱密度

设非平稳随机过程X(t)的自相关函数R_X（t₁,t₂)是两个时间点的函数，则$R_X(t_1,t_2)$的二维傅里叶变换
$$G_X({\omega }_1,{\omega }_2)={\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{R}_X(t_1,t_2)e^{-j({\omega }_2{t}_1-{\omega }_2{t}_2)}d{t}_{1}d{t}_2$$
为X(t)的广义功率谱密度，逆变换为
$$R_X(t_1,t_2)=\frac{1}{4{\pi }^2}{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }{G}_X({\omega }_1,{\omega }_2)e^{j({\omega }_1{t}_1-{\omega }_2{t}_2)}d{\omega }_{1}d{\omega }_2$$
广义谱的定义在数学上与平稳随机过程是相通的。假设$G_X({\omega }_1,{\omega }_2)$仅在两者相等处有值，两者不等时的结果是0，则
$$G_X({\omega }_1,{\omega }_2)=G_X({\omega }_1)\delta ({\omega }_1-{\omega }_2)$$这就是平稳情况下的功率谱。

时变功率谱

假设t₁=t,t₂=t+τ，则非平稳过程的自相关函数
$$R_X(t+\tau ,t)=R_X(t,\tau )=E[X(t+\tau )X(t)]$$
是t和$\tau$的函数，时变谱定义成
$$G_X(t,\omega )={\int }_{-\infty }^{\infty }{R}_X(t,\tau )e^{-j\omega \tau }d\tau$$
如果采用对称相关函数
$$R_X(t,\tau )=E[X(t+\tau /2)X(t-\tau /2)]$$
则这种时变谱称为维格纳-威利谱（WV谱）其中，$G_X(t,\omega )={\int }_{-\infty }^{\infty }{R}_X(t,\tau )e^{-j\omega \tau }d\tau$还可以写成
$$G_X(t,\omega )=E[W_X(t,\omega )]$$
其中
$$W_X(t,\omega )={\int }_{-\infty }^{\infty }[X(t+\tau /2)X(t-\tau /2)]e^{-j\omega \tau }d\tau$$
$W_X(t,\omega )$称为维格纳分布，非平稳随机过程X(t)的WV时变谱是该过程维格纳分布的数学期望。

例子

设噪声调制的振荡信号
$$X(t)=N(t)cos{\omega }_{0}t$$其中N(t)是平稳噪声，求X(t)的功率谱。
观察这个信号，发现不同时刻的信号（可能的）均值受到了三角函数的调制。X(t)是非平稳随机过程。求功率谱的方法一般就是求相关函数的傅里叶变换，而对于非平稳过程，有广义功率谱密度、时变功率谱等表示形式，但最常用的还是对即时相关函数做时间平均，然后进行傅里叶变换，如下。
$$R_X(\tau ,t)=\frac{1}{2}{R}_N(\tau )[cos{\omega }_0(2t+\tau )+cos{\omega }_0\tau ]$$
所谓的对时间函数求平均，指的是
$$\overline{R_X(\tau )}=\underset{T\to \infty }{\lim }\frac{1}{2T}{\int }_{-T}^T\frac{1}{2}{R}_N(\tau )[cos{\omega }_0(2t+\tau )+cos{\omega }_0\tau ]dt$$$$=\frac{1}{2}{R}_N(\tau )cos{\omega }_0\tau$$
对上述自相关函数求时间平均就得到了结果：
$$G_X(\omega )=\frac{1}{4}[G_N(\omega +{\omega }_0)+G_N(\omega -{\omega }_0)]$$

matlab互相关函数估计

参见文档https://www.mathworks.com/help/matlab/ref/xcorr.html

>> a=0.8;
>> sigma=2;
>> N=500;
>> u=randn(N,1);
>> x(1)=sigma*u(1)/sqrt(1-a^2);
>> for i=2:N
x(i)=a*x(i-1)+sigma*u(i);
end
>> R=xcorr(x,'coeff');
>> plot(R,'r','linewidth',2);

matlab功率谱估计

参照函数https://www.mathworks.com/help/signal/ref/periodogram.html

例如，估计两个正弦信号加正态白噪声的功率谱，信号
$$s(t)=cos2\pi f_{1}t+cos2\pi f_2(t)$$其中，$f_1=300Hz,f_2=310Hz$

>> t=0:0.001:0.3;
x=cos(2*pi*t*300)+cos(2*pi*t*310)+randn(size(t));
subplot(2,2,1);
periodogram(x,[],512,10000);
axis([0 500 -50 0]);
subplot(2,2,2);
window = hann(301);
periodogram(x,window,512,1000);
axis([0 500 -50 0]);
>> R=xcorr(x)/150000;
>> Pw=fft(R);
>> subplot(2,2,3);
>> f=(0:length(Pw)-1)*1000/length(Pw);
>> plot(f,10*log10(abs(Pw));
>>> plot(f,10*log10(abs(Pw)));;
>> axis([0 100 -50 0]);
>> grid on
>> subplot(2,2,4);
>> pwelvh(x,128,64,[],1000);
>>> pwelch(x,128,64,[],1000);
>> axis([0 500 -50 0]);

不加窗的周期图法如下图。就是先对离散序列做傅里叶变换，再将功率谱估计成频域能量的时间平均。
$$X(\omega )=\sum _{n=0}^{N-1}x(n)e^{-jn\omega }$$
则功率谱估计成是
$${\hat{G}}_X(\omega )=\frac{1}{N}|X(\omega ){|}^2$$

实际情况中数据截取的长度是有限的，数据的截断会导致一定截断误差的产生。减少截断效应，常采用数据加窗。下图是加汉宁窗后的周期图功率谱估计。

也可以采用传统的自相关函数法估计（先求出相关函数的估计，再对相关函数做傅里叶变换）

最后给出matlab另一个自带的韦尔奇功率谱估计。

例子：数字图像直方图

随机过程可以随时间变化，也可以看成随空间变化。数字图像可以看作随位置变化的随机序列。
数字图像灰度级的直方图，指的是反映图像中灰度级与出现这种灰度之间概率的图形。设R代表图像中像素灰度级，R的取值范围是[0,L-1]，L为总的灰度级数，
$$h(r_k)=n_{k}k=0,1,..,L-1$$
$r_k$：第k个灰度级，越往下越黑。n_k是具有灰度级$r_k$的像素度。直接想到可以对灰度做归一化，$$f(r_k)=\frac{n_k}{n}(k=0,1,...L-1)$$
f(r_k)是灰度级出现的概率估计。归一后就有了$$\sum _{k=0}^{L-1}f(r_k)=1$$
直方图可用于图像压缩、图像增强等处理技术中。下面是最最简单的处理：

>> I=imread('微信图片_20200427230532.jpg');
>> figure
subplot(1,2,1)
imshow(I)
subplot(1,2,2)
imhist(I,64)

调亮：I=I+30;直方图灰度级集中向高端移动。

均衡（从上面图的像素灰度级分布转换成服从均匀分布的灰度级。）
不妨先假设R是连续的，它（灰度级）归一化了，概率密度$f_R(r)$，对R做变换得到新变量S
$$S=T(R)$$
转移函数是0到1之间的单调函数。S的概率密度
$$f_s(s)=f_R(r)|\frac{dr}{ds}|$$
比较经典的是积分一下，相当于求分布函数
$$S=T(R)={\int }_0^R{f}_R(r)dr=F_R(R)$$
可以证明S在0到1服从均匀分布（利用随机变量的分布函数求解其反函数可得到任意分布随机数，随机变量的函数变换可获得任意分布随机数），$f_s(s)=1$。
这里把R变成离散的之后，就是把上面概率密度变成概率，把积分改成求和。R=r_k的概率是
P { R = r k } = f R ( r k ) = n k n P\{R=r_k\}=f_R(r_k)=\frac{n_k}{n} P{R=rk​}=fR​(rk​)=nnk​​
那同理，$$s_k=T(r_k)=\sum _{j=0}^k{f}_R(r_j)=\sum _{j=0}^k\frac{n_j}{n}$$
这样，图像像素由R变成S，就服从均匀分布了。matlab仍然已经写好了这个过程。

>> J=histeq(I);
>> figure
subplot(1,2,1)
imshow(J)
subplot(1,2,2)
imhist(J,64)