互功率谱

类似从自相关函数推出互相关函数，也可以从单个随机过程的功率谱推出互功率谱。两个随机过程X(t)和Y(t)，它们的互功率谱定义是

$G_{XY}(omega)=E[lim_{Trightarrow infty}frac{1}{2T}X_T(omega)Y_T^*(omega)]$

GXY​(ω)=E[T→∞lim​2T1​XT​(ω)YT∗​(ω)]

$X_T(omega)=int_{-T}^TX(t)e^{-jomega t}dt$

XT​(ω)=∫−TT​X(t)e−jωtdt

$Y_T(omega)=int_{-T}^TY(t)e^{-jomega t}dt$

YT​(ω)=∫−TT​Y(t)e−jωtdt
如果X(t)和Y(t)是联合平稳的，则
$R_{XY}(tau)$

RXY​(τ)绝对可积。则它们的互功率谱和互相关函数是一对傅里叶变换对：

$G_{XY}(omega)=int_{-infty}^{infty}R_{XY}(tau)e^{-jomega tau}dtau$

GXY​(ω)=∫−∞infty​RXY​(τ)e−jωτdτ

$R_{XY}(tau)=frac{1}{2pi}int_{-infty}^{infty}G_{XY}(omega)e^{jomega tau}domega$

RXY​(τ)=2π1​∫−∞∞​GXY​(ω)ejωτdω
互功率谱的性质：

$G_{XY}(omega)=G_{YX}^*(omega)$

GXY​(ω)=GYX∗​(ω)

$|G_{XY}(omega)|^2le G_X(omega)G_Y(omega)$

∣GXY​(ω)∣2≤GX​(ω)GY​(ω)

非平稳过程的功率谱

广义功率谱密度

设非平稳随机过程X(t)的自相关函数R_X（t₁,t₂)是两个时间点的函数，则
$R_X(t_1,t_2)$

RX​(t1​,t2​)的二维傅里叶变换

$G_X(omega_1,omega_2)=int_{-infty}^{+infty}int_{-infty}^{+infty}R_X(t_1,t_2)e^{-j(omega_2t_1-omega_2t_2)}dt_1dt_2$

GX​(ω1​,ω2​)=∫−∞+∞​∫−∞+∞​RX​(t1​,t2​)e−j(ω2​t1​−ω2​t2​)dt1​dt2​
为X(t)的广义功率谱密度，逆变换为

$R_X(t_1,t_2)=frac{1}{4pi^2}int_{-infty}^{+infty}int_{-infty}^{+infty}G_X(omega_1,omega_2)e^{j(omega_1t_1-omega_2t_2)}domega_1domega_2$

RX​(t1​,t2​)=4π21​∫−∞+∞​∫−∞+∞​GX​(ω1​,ω2​)ej(ω1​t1​−ω2​t2​)dω1​dω2​
广义谱的定义在数学上与平稳随机过程是相通的。假设
$G_X(omega_1,omega_2)$

GX​(ω1​,ω2​)仅在两者相等处有值，两者不等时的结果是0，则

$G_X(omega_1,omega_2)=G_X(omega_1)delta (omega_1-omega_2)$

GX​(ω1​,ω2​)=GX​(ω1​)δ(ω1​−ω2​)这就是平稳情况下的功率谱。

时变功率谱

假设t₁=t,t₂=t+τ，则非平稳过程的自相关函数

$R_X(t+tau,t)=R_X(t,tau)=E[X(t+tau)X(t)]$

RX​(t+τ,t)=RX​(t,τ)=E[X(t+τ)X(t)]
是t和
$tau$

τ的函数，时变谱定义成

$G_X(t,omega)=int_{-infty}^{infty}R_X(t,tau)e^{-jomega tau}dtau$

GX​(t,ω)=∫−∞∞​RX​(t,τ)e−jωτdτ
如果采用对称相关函数

$R_X(t,tau)=E[X(t+tau/2)X(t-tau/2)]$

RX​(t,τ)=E[X(t+τ/2)X(t−τ/2)]
则这种时变谱称为维格纳-威利谱（WV谱）其中，
$G_X(t,omega)=int_{-infty}^{infty}R_X(t,tau)e^{-jomega tau}dtau$

GX​(t,ω)=∫−∞∞​RX​(t,τ)e−jωτdτ还可以写成

$G_X(t,omega)=E[W_X(t,omega)]$

GX​(t,ω)=E[WX​(t,ω)]
其中

$W_X(t,omega)=int_{-infty}^{infty}[X(t+tau/2)X(t-tau/2)]e^{-jomega tau}dtau$

WX​(t,ω)=∫−∞∞​[X(t+τ/2)X(t−τ/2)]e−jωτdτ

$W_X(t,omega)$

WX​(t,ω)称为维格纳分布，非平稳随机过程X(t)的WV时变谱是该过程维格纳分布的数学期望。

例子

设噪声调制的振荡信号

$X(t)=N(t)cosomega_0t$

X(t)=N(t)cosω0​t其中N(t)是平稳噪声，求X(t)的功率谱。
观察这个信号，发现不同时刻的信号（可能的）均值受到了三角函数的调制。X(t)是非平稳随机过程。求功率谱的方法一般就是求相关函数的傅里叶变换，而对于非平稳过程，有广义功率谱密度、时变功率谱等表示形式，但最常用的还是对即时相关函数做时间平均，然后进行傅里叶变换，如下。

$R_X(tau,t)=frac{1}{2}R_N(tau)[cosomega_0(2t+tau)+cosomega_0tau]$

RX​(τ,t)=21​RN​(τ)[cosω0​(2t+τ)+cosω0​τ]
所谓的对时间函数求平均，指的是

$overline {R_X(tau)}=lim_{Trightarrowinfty}frac{1}{2T}int_{-T}^{T}frac{1}{2}R_N(tau)[cosomega_0(2t+tau)+cosomega_0tau]dt$

RX​(τ)​=T→∞lim​2T1​∫−TT​21​RN​(τ)[cosω0​(2t+τ)+cosω0​τ]dt
$=frac{1}{2}R_N(tau)cosomega_0tau$

=21​RN​(τ)cosω0​τ
对上述自相关函数求时间平均就得到了结果：

$G_X(omega)=frac{1}{4}[G_N(omega+omega_0)+G_N(omega-omega_0)]$

GX​(ω)=41​[GN​(ω+ω0​)+GN​(ω−ω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)=cos2pi f_1t+cos2pi f_2(t)$

s(t)=cos2πf1​t+cos2πf2​(t)其中，
$f_1=300Hz,f_2=310Hz$

f1​=300Hz,f2​=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^{-jnomega}$

X(ω)=n=0∑N−1​x(n)e−jnω
则功率谱估计成是

$hat G_X(omega)=frac{1}{N}|X(omega)|^2$

G^X​(ω)=N1​∣X(ω)∣2

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

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

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

例子：数字图像直方图

随机过程可以随时间变化，也可以看成随空间变化。数字图像可以看作随位置变化的随机序列。
数字图像灰度级的直方图，指的是反映图像中灰度级与出现这种灰度之间概率的图形。设R代表图像中像素灰度级，R的取值范围是[0,L-1]，L为总的灰度级数，

$h(r_k)=n_k~~k=0,1,..,L-1$

h(rk​)=nk​  k=0,1,..,L−1

$r_k$

rk​：第k个灰度级，越往下越黑。n_k是具有灰度级
$r_k$

rk​的像素度。直接想到可以对灰度做归一化，
$f(r_k)=frac{n_k}{n}~~(k=0,1,…L-1)$

f(rk​)=nnk​​  (k=0,1,...L−1)
f(r_k)是灰度级出现的概率估计。归一后就有了
$sum_{k=0}^{L-1}f(r_k)=1$

k=0∑L−1​f(rk​)=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)$

fR​(r)，对R做变换得到新变量S

$S=T(R)$

S=T(R)
转移函数是0到1之间的单调函数。S的概率密度

$f_s(s)=f_R(r)|frac{dr}{ds}|$

fs​(s)=fR​(r)∣dsdr​∣
比较经典的是积分一下，相当于求分布函数

$S=T(R)=int_0^Rf_R(r)dr=F_R(R)$

S=T(R)=∫0R​fR​(r)dr=FR​(R)
可以证明S在0到1服从均匀分布（利用随机变量的分布函数求解其反函数可得到任意分布随机数，随机变量的函数变换可获得任意分布随机数），
$f_s(s)=1$

fs​(s)=1。
这里把R变成离散的之后，就是把上面概率密度变成概率，把积分改成求和。R=r_k的概率是

$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}^kfrac{n_j}{n}$

sk​=T(rk​)=j=0∑k​fR​(rj​)=j=0∑k​nnj​​
这样，图像像素由R变成S，就服从均匀分布了。matlab仍然已经写好了这个过程。

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

随机过程应用：判断说话中的清音浊音

衡量标准综合如下：
语音片段的能量：短时能量通过大约长度为20ms的片段进行衡量。显然，浊音片段往往比清音片段有更高的能量，这从波形中就很容易看出。
过零率（zero crossing rate. ZCR）：
在昨天的例子中实际上已经实现了：

#在导入一段音频后，取一小段分析 n0= 57000 n1 = 57100 plt.figure(figsize = (14,5)) plt.plot(x[n0:n1]) plt.grid() 

计算过零点个数：

zero_crossings = librosa.zero_crossings(x[n0:n1],pad = False) print(sum(zero_crossings)) 

往往，浊音部分的过零率小 ，清音部分的过零率大（频率往往比较高）。而沉静片段的过零率往往和说话的过零率在一个量级，都比清音的低。当然，沉静片段的能量会小很多，所以可以区别开来。这里举一个例子，我说“东风来了，春天的脚步近了”，然后让它判断清音和浊音大概怎么分类。

首先读入文件并画图。读取前8000个样本点。

>> a=fread(fp,8000); >> plot(a) >> fp=fopen('春天的脚步近了.m4a'); fseek(fp,60244,-1); >> a=fread(fp,8000); >> plot(a) 

像刚才在python中做的一样，计算零交叉点的数目。8000确实有点多了，这里不妨取第5000到第7000的样本点进行处理。(已事先在音频软件中观测到5000到7000的样本是在清音阶段)

figure; for i=5000:5199 if (a(i)>128) && (a(i+1)<128)          x=x+1; else x=x;     end end disp('number of zero crossing='); disp(x); 

再在音频软件中定位到春天an处，（浊音）作出这附近的曲线：

fseek(fp,65042,-1); b=fread(fp,2000); plot(b); 

用相同的方法，得到零交叉点数目，可见大概只有上个数目的不到一半。

再在音频软件中定位到脚步的j处（清音)，作出这附近的曲线：

fseek(fp,6724,-1); c=fread(fp,200); plot(c); 

可以看出，它的振幅变化幅度较浊音更大。

归一化自协方差

$C_1=frac{sum_{n-1}^{N}S(n)S(n-1)}{sqrt{sum_{n-1}^N S^(n)sum _{n-0}^{N-1}S^2(n)}}$

C1​=∑n−1N​S(n)∑n−0N−1​S2(n)​∑n−1N​S(n)S(n−1)​
这个量规定了相邻音乐/语音片段的相关性。讲人话，C₁就可以理解成这时求的相邻单位时间延迟的自协方差。对于浊音，有大量低频信号能量，在浊音附近是的信号高度自相关的，这个值很大(讲人话，低频信号说明信号的乱起八糟的高频谐波占比比较少，信号比较“单纯”，相互的依赖性比较强)；而在清音区域自协方差就很小，在寂静区域会更小（讲人话，就是噪声基本上无法预测、是独立产生的）。与之类似，在语音信号中通常再引入一个叫谱倾斜的东西，定义成

$frac{sum_{i=1}^{N}s(i)s(i-1)}{sum_{i=1}s^2(i)}$

∑i=1​s2(i)∑i=1N​s(i)s(i−1)​
可以用matlab对“春天来了”做差分，求出它的谱倾斜度（spectrum tilt）（当然，需要事先录制”春天来了”4字）

>>  fp=fopen('录音 (5).m4a'); %fseek(fp,800,-1); a=fread(fp,400,'short'); a=a-128; subplot(2,1,1);plot(a); xlabel('sample no.');ylabel('amplitude'); for j=1:200, fseek(fp,4*j,-1);   a=fread(fp,100);   a=a-128; sum(j)=0; for i=2:100, sum(j)=sum(j)+(a(i)*a(i-1));   end  sum1(j)=0; for i=2:100, if a(i)==0, a(i)=0.1;    end   sum1(j)=sum1(j)+a(i)*a(i);   end  s(j)=sum(j)/sum1(j); end subplot(2,1,2);plot(s); 

得到的是下图。

可以看出，由于我说的比较波澜不惊，波形很规整的了。音频信号的自相关系数（谱倾斜度和自相关函数实际上差不多是同一种量度）在浊音处是十分高的；由于春和天是连在一起的，所以得到的相关性都比较大（接近1）；而后来说来字的时候，音调和共振峰都发生了比较强烈的改变，所以它的相关性急剧下降。