互功率谱
类似从自相关函数推出互相关函数,也可以从单个随机过程的功率谱推出互功率谱。两个随机过程X(t)和Y(t),它们的互功率谱定义是
G
X
Y
(
ω
)
=
E
[
lim
T
→
∞
1
2
T
X
T
(
ω
)
Y
T
∗
(
ω
)
]
G_{XY}(\omega)=E[\lim_{T\rightarrow \infty}\frac{1}{2T}X_T(\omega)Y_T^*(\omega)]
GXY(ω)=E[T→∞lim2T1XT(ω)YT∗(ω)]
X
T
(
ω
)
=
∫
−
T
T
X
(
t
)
e
−
j
ω
t
d
t
X_T(\omega)=\int_{-T}^TX(t)e^{-j\omega t}dt
XT(ω)=∫−TTX(t)e−jωtdt
Y
T
(
ω
)
=
∫
−
T
T
Y
(
t
)
e
−
j
ω
t
d
t
Y_T(\omega)=\int_{-T}^TY(t)e^{-j\omega t}dt
YT(ω)=∫−TTY(t)e−jωtdt
如果X(t)和Y(t)是联合平稳的,则
R
X
Y
(
τ
)
R_{XY}(\tau)
RXY(τ)绝对可积。则它们的互功率谱和互相关函数是一对傅里叶变换对:
G
X
Y
(
ω
)
=
∫
−
∞
i
n
f
t
y
R
X
Y
(
τ
)
e
−
j
ω
τ
d
τ
G_{XY}(\omega)=\int_{-\infty}^{infty}R_{XY}(\tau)e^{-j\omega \tau}d\tau
GXY(ω)=∫−∞inftyRXY(τ)e−jωτdτ
R
X
Y
(
τ
)
=
1
2
π
∫
−
∞
∞
G
X
Y
(
ω
)
e
j
ω
τ
d
ω
R_{XY}(\tau)=\frac{1}{2\pi}\int_{-\infty}^{\infty}G_{XY}(\omega)e^{j\omega \tau}d\omega
RXY(τ)=2π1∫−∞∞GXY(ω)ejωτdω
互功率谱的性质:
G
X
Y
(
ω
)
=
G
Y
X
∗
(
ω
)
G_{XY}(\omega)=G_{YX}^*(\omega)
GXY(ω)=GYX∗(ω)
∣
G
X
Y
(
ω
)
∣
2
≤
G
X
(
ω
)
G
Y
(
ω
)
|G_{XY}(\omega)|^2\le G_X(\omega)G_Y(\omega)
∣GXY(ω)∣2≤GX(ω)GY(ω)
非平稳过程的功率谱
广义功率谱密度
设非平稳随机过程X(t)的自相关函数RX(t1,t2)是两个时间点的函数,则
R
X
(
t
1
,
t
2
)
R_X(t_1,t_2)
RX(t1,t2)的二维傅里叶变换
G
X
(
ω
1
,
ω
2
)
=
∫
−
∞
+
∞
∫
−
∞
+
∞
R
X
(
t
1
,
t
2
)
e
−
j
(
ω
2
t
1
−
ω
2
t
2
)
d
t
1
d
t
2
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(ω2t1−ω2t2)dt1dt2
为X(t)的广义功率谱密度,逆变换为
R
X
(
t
1
,
t
2
)
=
1
4
π
2
∫
−
∞
+
∞
∫
−
∞
+
∞
G
X
(
ω
1
,
ω
2
)
e
j
(
ω
1
t
1
−
ω
2
t
2
)
d
ω
1
d
ω
2
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_1t_1-\omega_2t_2)}d\omega_1d\omega_2
RX(t1,t2)=4π21∫−∞+∞∫−∞+∞GX(ω1,ω2)ej(ω1t1−ω2t2)dω1dω2
广义谱的定义在数学上与平稳随机过程是相通的。假设
G
X
(
ω
1
,
ω
2
)
G_X(\omega_1,\omega_2)
GX(ω1,ω2)仅在两者相等处有值,两者不等时的结果是0,则
G
X
(
ω
1
,
ω
2
)
=
G
X
(
ω
1
)
δ
(
ω
1
−
ω
2
)
G_X(\omega_1,\omega_2)=G_X(\omega_1)\delta (\omega_1-\omega_2)
GX(ω1,ω2)=GX(ω1)δ(ω1−ω2)这就是平稳情况下的功率谱。
时变功率谱
假设t1=t,t2=t+τ,则非平稳过程的自相关函数
R
X
(
t
+
τ
,
t
)
=
R
X
(
t
,
τ
)
=
E
[
X
(
t
+
τ
)
X
(
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
,
ω
)
=
∫
−
∞
∞
R
X
(
t
,
τ
)
e
−
j
ω
τ
d
τ
G_X(t,\omega)=\int_{-\infty}^{\infty}R_X(t,\tau)e^{-j\omega \tau}d\tau
GX(t,ω)=∫−∞∞RX(t,τ)e−jωτdτ
如果采用对称相关函数
R
X
(
t
,
τ
)
=
E
[
X
(
t
+
τ
/
2
)
X
(
t
−
τ
/
2
)
]
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
,
ω
)
=
∫
−
∞
∞
R
X
(
t
,
τ
)
e
−
j
ω
τ
d
τ
G_X(t,\omega)=\int_{-\infty}^{\infty}R_X(t,\tau)e^{-j\omega \tau}d\tau
GX(t,ω)=∫−∞∞RX(t,τ)e−jωτdτ还可以写成
G
X
(
t
,
ω
)
=
E
[
W
X
(
t
,
ω
)
]
G_X(t,\omega)=E[W_X(t,\omega)]
GX(t,ω)=E[WX(t,ω)]
其中
W
X
(
t
,
ω
)
=
∫
−
∞
∞
[
X
(
t
+
τ
/
2
)
X
(
t
−
τ
/
2
)
]
e
−
j
ω
τ
d
τ
W_X(t,\omega)=\int_{-\infty}^{\infty}[X(t+\tau/2)X(t-\tau/2)]e^{-j\omega \tau}d\tau
WX(t,ω)=∫−∞∞[X(t+τ/2)X(t−τ/2)]e−jωτdτ
W
X
(
t
,
ω
)
W_X(t,\omega)
WX(t,ω)称为维格纳分布,非平稳随机过程X(t)的WV时变谱是该过程维格纳分布的数学期望。
例子
设噪声调制的振荡信号
X
(
t
)
=
N
(
t
)
c
o
s
ω
0
t
X(t)=N(t)cos\omega_0t
X(t)=N(t)cosω0t其中N(t)是平稳噪声,求X(t)的功率谱。
观察这个信号,发现不同时刻的信号(可能的)均值受到了三角函数的调制。X(t)是非平稳随机过程。求功率谱的方法一般就是求相关函数的傅里叶变换,而对于非平稳过程,有广义功率谱密度、时变功率谱等表示形式,但最常用的还是对即时相关函数做时间平均,然后进行傅里叶变换,如下。
R
X
(
τ
,
t
)
=
1
2
R
N
(
τ
)
[
c
o
s
ω
0
(
2
t
+
τ
)
+
c
o
s
ω
0
τ
]
R_X(\tau,t)=\frac{1}{2}R_N(\tau)[cos\omega_0(2t+\tau)+cos\omega_0\tau]
RX(τ,t)=21RN(τ)[cosω0(2t+τ)+cosω0τ]
所谓的对时间函数求平均,指的是
R
X
(
τ
)
‾
=
lim
T
→
∞
1
2
T
∫
−
T
T
1
2
R
N
(
τ
)
[
c
o
s
ω
0
(
2
t
+
τ
)
+
c
o
s
ω
0
τ
]
d
t
\overline {R_X(\tau)}=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}\frac{1}{2}R_N(\tau)[cos\omega_0(2t+\tau)+cos\omega_0\tau]dt
RX(τ)=T→∞lim2T1∫−TT21RN(τ)[cosω0(2t+τ)+cosω0τ]dt
=
1
2
R
N
(
τ
)
c
o
s
ω
0
τ
=\frac{1}{2}R_N(\tau)cos\omega_0\tau
=21RN(τ)cosω0τ
对上述自相关函数求时间平均就得到了结果:
G
X
(
ω
)
=
1
4
[
G
N
(
ω
+
ω
0
)
+
G
N
(
ω
−
ω
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
)
=
c
o
s
2
π
f
1
t
+
c
o
s
2
π
f
2
(
t
)
s(t)=cos2\pi f_1t+cos2\pi f_2(t)
s(t)=cos2πf1t+cos2πf2(t)其中,
f
1
=
300
H
z
,
f
2
=
310
H
z
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
(
ω
)
=
∑
n
=
0
N
−
1
x
(
n
)
e
−
j
n
ω
X(\omega)=\sum_{n=0}^{N-1}x(n)e^{-jn\omega}
X(ω)=n=0∑N−1x(n)e−jnω
则功率谱估计成是
G
^
X
(
ω
)
=
1
N
∣
X
(
ω
)
∣
2
\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(r_k)=n_k~~k=0,1,..,L-1
h(rk)=nk k=0,1,..,L−1
r
k
r_k
rk:第k个灰度级,越往下越黑。nk是具有灰度级
r
k
r_k
rk的像素度。直接想到可以对灰度做归一化,
f
(
r
k
)
=
n
k
n
(
k
=
0
,
1
,
.
.
.
L
−
1
)
f(r_k)=\frac{n_k}{n}~~(k=0,1,...L-1)
f(rk)=nnk (k=0,1,...L−1)
f(rk)是灰度级出现的概率估计。归一后就有了
∑
k
=
0
L
−
1
f
(
r
k
)
=
1
\sum_{k=0}^{L-1}f(r_k)=1
k=0∑L−1f(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
)
f_R(r)
fR(r),对R做变换得到新变量S
S
=
T
(
R
)
S=T(R)
S=T(R)
转移函数是0到1之间的单调函数。S的概率密度
f
s
(
s
)
=
f
R
(
r
)
∣
d
r
d
s
∣
f_s(s)=f_R(r)|\frac{dr}{ds}|
fs(s)=fR(r)∣dsdr∣
比较经典的是积分一下,相当于求分布函数
S
=
T
(
R
)
=
∫
0
R
f
R
(
r
)
d
r
=
F
R
(
R
)
S=T(R)=\int_0^Rf_R(r)dr=F_R(R)
S=T(R)=∫0RfR(r)dr=FR(R)
可以证明S在0到1服从均匀分布(利用随机变量的分布函数求解其反函数可得到任意分布随机数,随机变量的函数变换可获得任意分布随机数),
f
s
(
s
)
=
1
f_s(s)=1
fs(s)=1。
这里把R变成离散的之后,就是把上面概率密度变成概率,把积分改成求和。R=rk的概率是
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
)
=
∑
j
=
0
k
f
R
(
r
j
)
=
∑
j
=
0
k
n
j
n
s_k=T(r_k)=\sum_{j=0}^{k}f_R(r_j)=\sum_{j=0}^k\frac{n_j}{n}
sk=T(rk)=j=0∑kfR(rj)=j=0∑knnj
这样,图像像素由R变成S,就服从均匀分布了。matlab仍然已经写好了这个过程。
>> J=histeq(I);
>> figure
subplot(1,2,1)
imshow(J)
subplot(1,2,2)
imhist(J,64)