循环平稳过程以及信号处理理论

循环平稳过程以及信号处理理论
绪论
• 通信、遥测、雷达、声呐等系统中许多信号，其统计特征参数是时间变化的，这类信号称为循环平稳信号(cyclostationary signal)
• 例如调制信号，雷达扫描信号，还有一些自然的，如水文数据，海洋数据，人体心电图等都具有循环平稳性质。

• W. A. Gardner*的谱相关理论是标志循环平稳信号处理理论的成熟，其数学工具是循环相关函数和循环谱相关函数。

• *W. A. Gardner, L. E. Franks, Characterization of cyclostationary random signal processes, IEEE Trans Information Theory, 21: 4-14, 1975.
• F. Chapeau-Blondeau, X. Godivier; "Theory of stochastic resonance in signal
transmission by static nonlinear systems"; Physical Review E 55, 1478-1495 (1997).
• X. Godivier, F. Chapeau-Blondeau; "Noise-assisted signal transmission by a nonlinear
electronic comparator: Experiment and Theory"; Signal Processing 56, 293-303 (1997).
• F. Duan, F. Chapeau-Blondeau, D. Abbott; "Noise-enhanced SNR gain in parallel array of
bistable oscillators"; Electronics Letters 42, 1008-1009 (2006).
2.1一般理论框架（动态静态系统都适合）
强调我们的系统划分规则静态指无记忆系统，而动态指有记忆的系统。

这里设任意一系统的输入为)()()(t t s t x η+=，)(t s 表示周期为s T 的周期信号，而)(t η是稳态随机噪声。

我们把系统输出)(t y 看成是它的非稳态均值)]([t y E 与围绕均值的稳定波动
)(ˆt y
的和，即 )]([)(ˆ)(t y E t y
t y += 2-1 由于输入信号)(t s 的周期性，系统输出)(t y 一般也是周期为s T 的循环平稳信号，非稳态均值)]([t y E 是周期s T 的确定性信号，那么引入傅里叶变换系数n Y
这里整数,...2,1,0=n ，表示了谐波频率倍数。

为计算自相关函数，固定时间t 和时间延迟τ，得二阶非稳态相关函数
)]([)]([)](~)(~[)]()([τττ+++=+t y E t y E t y t y E t y t y E 2-3
此期望函数含有两个变量时间t 和时间延迟τ，同样具有周期性。

为构造稳态的自相关函数，我们进行时间平均可以得到)(t
y 的自相关函数)(τyy R
其中，自协方差)(τyy C 为
那么，依据维纳辛钦定理，)(t y 的功率谱密度)(v P yy ,即是)
(τyy R 的傅里叶变换
其中
τ
τ
τ
ττπτ
πππτπτππτπd e dte
e
Y e
Y T d e dt e
Y e
Y T d e dt t y E t y E T v j j k
t
j k T n
t
j n s
v j k
t j k T n
t
j n s
v j T s
s
k s
k s
s
n s
T k s
s
T n s
2220
22)
(20220
11)]([)]([1-∞
∞--+∞
∞--∞
∞
-∑⎰⎰∑∑⎰
⎰
∑⎰
⎰
==+
由Parseval 定理，当k n =时，
*20
2
2||||1
n n n T t j n s
Y Y Y dt e
Y T s
T n ==⎰
π
而当k n ≠时，t j T
n e
π2与t
j k
e
π2是正交基，
01
22=⎰
dt e
e
T s
T k
T n T t
j t
j s
ππ
而对于信号t j T
n e
π2的Fourier 变换为
)(22T n
v j j v d e
e
n
-=-∞
∞
-⎰
δττ
πτπ
由此得出2-6式。

由式2-6可以看出)(t y 的功率谱密度为在宽带噪声谱)]([τyy C F 背景下，在
s T /1的整数倍值上（谐频上）叠加了相关谱线。

其次分析)(τyy C ，由于自协方差
0)](~)(~[lim ||=+∞
→ττt y t y E
那么
由于)(t y 的非稳态方差)](ˆ)(ˆ[)](var[t y t y
E t y =，进行时间平均则为)](var[)](var[1
)0(0
t y dt t y T C T
s
yy ==
⎰
，那么)(τyy C 可以写为
2-7
)(τh =)(τyy C /)0(yy C 为归一化稳态自相关函数，是偶函数，其Fourier 变换为)()]([v H h F =τ。

那么2-6式可以重新写为
输出信噪比)/(s out T n R 定义为谱线s T n /上的信号功率与以谱线s T n /为中心的附近频带B ∆内噪声功率的比值，即
2.2静态系统的离散实现
考虑无记忆的静态系统，其输入输出转换为
)],()([)(t t s g t y η+= 2-10
这里g 是任意实函数，设)(t η是白噪声，概率密度为ηf ，分布函数为
⎰
∞
-=u
dx x f u F )()(ηη。

噪声的自相关函数为
)(2)]()([)(τδτηητηηD t t E R =+=
这样导致噪声具有无限大的功率)0(ηηR ，实际上，白噪声近似有一个小但是非零的相关时间c τ，这样其功率很大但是有限，使得D R c 2~)0(τηη条件成立。

这里对于输入为白噪声来说，数值模拟中，若采样时间步长为s T t <<∆，
s T t N =∆，对于离散采样白噪声)(t j ∆η，其自相关函数变为
)()]()([)(2
t k t t k t j t j E t k R ∆∆=∆+∆∆=∆δσηηηηη 2-11
离散的Dirac 函数定义为
⎩⎨
⎧≠=∆=∆0,
00
,/1)(k k t t k δ 那么功率谱密度为
t D ∆=2
2ησ
由于噪声为白噪声，那么静态非线性)(x g 的输出)(t y 和)(τ+t y 也具有非相关性，0≠τ。

那么，当0≠k 时二阶矩
)]([)]([)]()([t k t j y E t j y E t k t j y t j y E ∆+∆∆=∆+∆∆
当0=k 时，二阶矩
)]([)](~[)]()([22t j y E t j y E t j y t j y E ∆+∆=∆∆
在任意给定时刻t j t ∆=，输入)()(t t s η+的概率密度函数可以表示为
)(s u f -η，那么我们可以求解一阶矩
⎰∞
∞--=du t s u f u g t y E )]([)()]([η
二阶矩
⎰∞
∞
--=du t s u f u g t y E )]([)()]([22η
非稳态方差
)]([)]([)](ˆ)(ˆ[)](var[22t y E t y E t y t y
E t y -==， 因此，我们将)]()([t k t j y t j y E ∆+∆∆可以写成
)]([)]([)()](var[)]()([t k t j y E t j y E t k t t j y t k t j y t j y E ∆+∆∆+∆∆∆=∆+∆∆δ 2-12
对于时间t j t ∆=求时间平均得到稳态自相关函数
)]
([)]([1)()var()]
([)]([1
)()](var[1
)(1
1
1
1
1
1
t k t j y E t j y E N
t k t y t k t j y E t j y E N
t k t t j y N
t k R N j N j N j y ∆+∆∆+
∆∆=∆+∆∆+
∆∆∆=
∆∑
∑
∑
-=-=-=δδ2-13
这里时间积分
∑∑⎰-=-==∆*∆=*1
1
110*1
11N j N j T N
t t N dt T ，
因此稳态自协方差
)()var()()var()(t k h y t k t y t k C y ∆=∆∆=∆δ
同样地，为了观测频域特征，将均值的Fourier 系数表示为
∑-=⎥⎦⎤⎢⎣⎡
-∆=
1
1
2exp )]([1N j n N jn i t j y E N
Y π 这里离散周期Fourier 变换
]2exp[*1]2exp[1]2exp[111110N
jn
i N t t j t N n i t N dt t T n i T N j N j T πππ-=∆∆∆-*∆=-*∑∑⎰-=-= 同样将自相关函数)(t k R y ∆进行离散Fourier 变换
t MN
kl
i t k R t k R DFT MN MN
k y y ∆-∆=
∆∑
--=1
]22exp[)()]([π
这里Fourier 变换为
t t k t
MN l
i d v i MN MN
k ∆∆∆-=
-∑
⎰
--=∞
∞
-1]22exp[*]2exp[*π
ττπ
这里无穷积分上下限用MT 2周期来代替，那么频率分辨率为
)2/(1)2/(1t MN MT v ∆==∆，]1,[--∈MN MN l ，v l v ∆=。

那么
v Y t y Y Y t MN t y T n P t MN jl
jl kl i t k t j y E t j y E N
t MN
kl
i t k y t MN kl
i t k R t k R DFT n n n MN MN
k N j MN MN k MN MN
k y y ∆+∆=∆+∆=⎪⎭
⎫
⎝⎛∆-+-∆+∆∆+
∆-∆=∆-∆=
∆∑
∑∑
∑
--=-=--=--=1||]var[2]var[]22exp[)]([)]([1]22exp[)(]var[]22exp[)()]([2
*11
1
1
ππδπ
2-14
第一项是由于只有0=k 时，t t k ∆=∆/1)(δ，后面的项计算参考文献[X
Godivier, F. Chapeau-Blondeau, Noise-assisted signal transmission in a nonlinear electronic comparator, Signal Processing, 56:293-303, 1997]中利用
∑--=⎥⎦⎤⎢⎣⎡
-=
1
22exp MN MN
k MN kl i W π的一些性质证明。

由公式2-14可以看出，功率谱为一个常数背景t y ∆]var[上，谱线T /1以及其谐频T n /。

那么信噪比定义为谱线T n /上信号的功率除以噪声在频带
B ∆中所包含的功率，即
实际运算中可取)/(1/1t N T B ∆==∆，那么N B t /1=∆∆，代入得到
这里一个周期内采样N 个，一般100,1000，采样点多些。

2.3 信噪比的数值实现
这里无论是有记忆的动态系统还是无记忆的静态系统，我们探讨一类计算机计算的数值方法。

1）再次回顾，把系统输出)(t y 看成是它的非稳态均值)]([t y E 与围绕均值的稳定
波动)(ˆt y
的和，即 )]([)(ˆ)(t y E t y
t y += 引入傅里叶变换系数n Y
我们只有求)]([t j y E ∆离散值，并且是一个周期内t T N ∆=/个，利用样本来求样
y 1 y 2 y 3 y N y N+1 ……
本均值)]([t j y E ∆可以转化为求取时间均值
N j ,2,1= 2-16
这里考虑到动态系统的输出稳定性，把)](t y ∆到)(t N y ∆第一个周期的输出给舍弃掉了。

对于])[(],)[(t j KN y t j N y ∆+∆+ 这N 个值进行时间平均，他们都是
)]([)(ˆ)(t y E t y
t y +=，平均化之后，0/)(ˆ→∑K t y ，而把)]([t j y E ∆浮现出来。

这样Fourier 系数近似计算为
∑-=⎥⎦⎤⎢⎣
⎡
-∆=112exp )(1N j n N jn i t j y N Y π 2-17
这里整数,...2,1,0=n ，表示了谐波频率倍数。

2）为计算自相关函数，固定时间t j t ∆=和时间延迟t i ∆=τ，得二阶非稳态相关函数)]()([t i t j y t j y E ∆+∆∆
同样利用时间平均方法求
])[(])[()()(1t i t j kN y t j kN y t i t j y t j y K
k ∆+∆+∆+=∆+∆∆∑= 2-18
自协方差)(τyy C 为
)]([)]([)]()([)](~)(~[τττ+-+=+t y E t y E t y t y E t y t y E
那么根据2-16和2-18，可得
[]
∑=∆+∆∆-∆+∆∆=∆N
j yy t i t j y t j y t i t j y t j y N t i C 1
)(*)()()(1)( 2-19
那么，依据维纳辛钦定理，)(t y 的功率谱密度
对于Fourier 系数n Y 2-17式已经求出，那么只要求出)]([t i C DFT y ∆即可，注意其是偶函数，且有周期性，一个周期内衰减很快，因此只计算一个周期内即可
t N i j t i C N t i C DFT N i y y ∆⎥⎦⎤⎢⎣
⎡
-∆=∆∑=π2exp )(1)]([1
那么输出信噪比即可计算出来。