第5章随机信号分析 前面几章介绍的信号都是确定性信号,本章主要介绍的是随机信号。随机信号和确定性信号不同,它不能通过一个确切的数学公式来描述,也不能被准确地预测。因此,我们只能在统计的意义上来研究随机信号。随机信号的分析和处理方法与确定性信号有较大的差异。 本章首先简单介绍随机信号的定义和分类,然后重点讨论随机信号的相关分析和功率谱估计,最后介绍随机信号通过线性移不变系统的响应。 5.1随机信号简介 反映随机物理现象的信号不能用精确的数学关系式来描述,因为对这种现象的每一次观测都是不一样的,即任意一次观测只代表许多可能产生的结果之一,这种数据就是随机信号。表示随机现象的单个时间历程,称为样本函数,在有限时间区间上观测时,称为样本记录。随机现象可能产生的全部样本函数的集合,称为随机过程。 随机过程可分为平稳过程和非平稳过程两类,平稳随机过程可进一步分为各态遍历和非各态遍历两类; 非平稳随机过程又可进一步分为一般非平稳随机过程和瞬变随机过程。 5.1.1平稳随机过程 若把某种物理现象看作一个随机过程,这种现象在任何时刻的特性就可以用随机过程样本函数集合的平均值来描述。例如,图5.1.1表示一随机过程的样本函数集合(也称总体),则随机过程在某一时刻t1的均值,就是将总体中t1时刻的各样本函数的瞬时值相加,然后除以样本函数的个数,这就是以后经常用到的总体平均概念。类似地,随机过程两个不同时刻之值的相关性也称为自相关函数,可以由t1和 t1+τ两时刻瞬时值乘积的总体平均得到[1]。 图5.1.1组成随机过程的样本函数总体 我们一般用符号{}表示样本函数的总体。于是,随机过程{x(t)}的均值μx(t1)和自相关函数Rx(t1,t1+τ)分别为 μx(t1)= limN→∞1N∑Nk=1xk(t1)(5.1.1) Rx(t1,t1+τ)= limN→∞1N∑Nk=1xk(t1)xk(t1+τ)(5.1.2) 其中,在最后求和时,假定各样本函数都有等可能性。 在一般情况下,μx(t1)和Rx(t1,t1+τ)都随t1改变而变化,此时随机过程{x(t)}为非平稳的。在特殊情况下,若μx(t1)和Rx(t1,t1+τ)不随t1改变而变化,则称随机过程{x(t)}为弱平稳的或广义平稳的。对于弱平稳过程,均值是常数,自相关函数仅与时间位移 τ有关,即 μx(t1)=μx (5.1.3) Rx(t1,t1+τ)=Rx(τ)(5.1.4) 5.1.2各态遍历随机过程 在大多数情况下,可以用总体中某样本函数的时间平均来确定平稳随机过程的特性。例如,图5.1.1所示的随机过程中的第 k 个样本函数,它的均值μx(k)和Rx(τ,k)分别为 μx(k)= limT→∞1T∫T0xk(t)dt (5.1.5) Rx(τ,k)= limT→∞1T∫T0xk(t)xk(t+τ)dt (5.1.6) 如果随机过程{x(t)}是平稳的,而且用不同样本函数计算式(5.1.5)和式(5.1.6)中μx(k)和Rx(τ,k)的结果都一样,则称此随机过程为各态遍历的。对于各态遍历的随机过程,按时间平均的均值和自相关函数,以及所有其他按时间平均的量都等于相应的随机过程总体平均值,即 μx(k)=μx,Rx(τ,k)=Rx(τ) 各态遍历随机过程的所有特性可以用单个样本函数上的时间平均来描述,因此,各态遍历随机过程显然是随机过程中很基本的一类。实际上,表示平稳物理现象的随机信号一般是近似各态历经的,在大多数情况下,可以用单个观测到的时间历程记录来测定平稳随机现象的总体特征。 5.1.3非平稳随机过程 非平稳随机过程包括所有不满足平稳性要求的随机过程。非平稳随机过程的特性一般是随时间而变化的,因而只能用组成过程的样本函数的总体瞬时平均来确定。在实际中,由于不容易得到足够数量的样本记录来精确地测量总体平均性质,因此,对非平稳随机过程的测试和分析是比较困难的。 5.2随机信号的相关分析 5.2.1随机信号的自相关函数及其应用 1. 定义 1) 一般随机信号 若 X(t)为连续随机信号,当t为固定值时,则其概率分布函数 PX(x,t)为 PX(x,t)=∫x-∞pX(x;t)dx (5.2.1) 概率密度 pX(x,t)为 pX(x;t)=PX(x,t)/x (5.2.2) 均值μX(t)为 μX(t)=E{X(t)}=∫∞-∞xpX(x;t)dx (5.2.3) μX(t)也称为随机过程的数学期望,它是一个随机过程各个实现的平均函数,随机过程就在它的附近起伏变化。均方值 D2X(t)为 D2X(t)=E{X2(t)}=∫∞-∞x2pX(x;t)dx (5.2.4) 式中,E{·}表示求均值运算。D2X(t)也称为随机变量 X(t)的二阶原点矩。方差σ2X(t)为 σ2X(t)=E{[X(t)-μX(t)]2}(5.2.5) σ2X(t)也可以表示为 σ2X(t)=E{X2(t)}-μ2X(t)(5.2.6) σ2X(t)也称为随机变量的二阶中心矩,它描述了随机过程中诸样本函数围绕数学期望μX(t)的分散程度。 定义自相关函数为 RX(t1,t2)=E{X(t1)X(t2)}=∫∞-∞∫∞-∞x1x2pX(x1,x2;t1,t2)dx1dx2 (5.2.7) 自相关函数就是用来描述随机过程 X(t)任意两个不同时刻状态之间相关性的重要数字特征,也是 X(t)在两个不同时刻的状态之间的混和矩,它反映了 X(t)在两个不同时刻的状态之间的统计关联程度。 当 t1=t2=t 时,RX(t,t)=E{X2(t)}。由式(5.2.6)可得 RX(t,t)=σ2X(t,t)+μ2X(t)(5.2.8) 若 X(n)为离散随机信号,其中每一次的实现为 x(n,i),i=1,2,…,N,N→∞,则X(n)的均值μX(n)为 μX(n)=E{X(n)}= limN→∞1N∑Ni=1x(n,i)(5.2.9) 方差σ2X(n)为 σ2X(n)=E{X(n)-μX(n)2}= limN→∞1N∑Ni=1x(n,i)-μX(n)2(5.2.10) 均方值D2X 为 D2X=E{X(n)2}= limN→∞1N∑Ni=1x(n,i)2(5.2.11) 定义 X(n)的自相关函数为 RX(n1,n2)=E{X*(n1)X(n2)}= limN→∞∑Ni=1x*(n1,i)x(n2,i) (5.2.12) 随机信号的自相关函数RX(n1,n2)描述了信号 X(n)在 n1,n2 这两个时刻的互相关系,它是一个重要的统计量。若 n1=n2=n,则 RX(n1,n2)=E{X(n)2}=D2X(n)(5.2.13) 2) 广义平稳随机信号 若 X(n)为广义平稳随机信号,其均值与时间 n 无关,自相关函数RX(n1,n2)和 n1,n2 的起点无关,而仅与 n1,n2 之差有关,即其均值为 μX(n)=μX=E{X(n)}(5.2.14) 自相关函数为 RX(n1,n2)=RX(m)=E{X*(n)X(n+m)},m=n2-n1 (5.2.15) 广义平稳随机信号是一类重要的随机信号。在实际中,大多数的随机信号都可以认为是广义平稳随机信号,这样,也可使问题得以简化[2]。 例5.2.1设随机相位正弦序列为 X(n)=Asin(2πfnTs+Φ) 式中,A 和 f 均为常数; Ts 是采样时间间隔; Φ是一个随机变量,在 0~2π 内服从均匀分布,即 p(φ)=12π,0≤φ≤2π 0,其他 显然,对应Φ的一个取值,可以得到一条正弦曲线,因为Φ在 0~2π 内的取值是随机的,所以,其每一个样本x(n)都是一个正弦信号。求其均值及其自相关函数,并判断其平稳性[3]。 解根据随机信号的定义,X(n)的均值和自相关函数分别是 μX(n)=E{Asin(2πfnTs+Φ)}=∫2π0Asin(2πfnTs+φ)12πdφ =0 RX(n1,n2)=E{A2sin(2πfn1Ts+Φ)sin(2πfn2Ts+Φ)} =A22π∫2π0sin(2πfn1Ts+φ)sin(2πfn2Ts+φ)dφ =A22cos[2πf(n2-n1)Ts] 由于 μX(n)=μX=0 及 RX(n1,n2)=RX(n2-n1)=RX(m)=A22cos(2πfmTs) 所以,随机相位正弦波是广义平稳的。 ■ 3) 各态遍历随机信号 随机信号的均值和自相关函数等都是建立在集合平均的意义上的,为了精确地求出自相关函数需要知道 x(n,i),i=1,2,…,∞,这在实际中是不现实的。我们在实际工作中往往只能得到 X(n)的一次实验记录,即一个样本函数,对于一部分平稳信号可以用一次实验记录来代替一族记录计算 X(n)的自相关函数。对一平稳信号 X(n),如果它的所有样本函数在某一固定时刻的一阶和二阶统计特性与单一样本函数在长时间内的统计特性是一致的,则称 X(n)为各态遍历信号,其意义是,单一样本函数随时间变化的过程可以包括该信号所有样本函数的取值经历。这样我们可以重新定义均值和自相关函数如下: μX=E{X(n)}= limN→∞12N+1∑Nn=-Nx(n)=μx(5.2.16) RX(m)=E{X(n)X(n+m)} = limN→∞12N+1∑Nn=-Nx(n)x(n+m)=Rx(m)(5.2.17) 上面两式右边都是使用的单一样本函数x(n)来求出μx 和Rx(m)的,因此,称为时间平均。对各态遍历信号,其一阶和二阶的集合平均等于相应的时间平均。 例5.2.2讨论例5.2.1随机相位正弦波的各态遍历性。 解对 X(n)=Asin(2πfnTs+Φ),其单一的时间样本为 x(n)=Asin(2πfnTs+φ) φ 为一常数,对 X(n)进行时间平均得 μx= limN→∞12N+1∑Nn=-NAsin(2πfnTs+φ)=0=μX Rx(m)= limN→∞12N+1∑Nn=-NA2sin(2πfnTs+φ)sin[2πf(n+m)Ts+φ] = limN→∞12N+1∑Nn=-NA22{cos(2πfmTs)-cos[2πf(n+n+m)Ts]} 由于上式是对 n 求和,故求和号中的第一项与 n 无关,而第二项应等于0,所以 Rx(m)=A22cos(2πfmTs)=RX(m) 这与例5.2.1按集合平均求出的结果一样,所以,随机相位正弦波既是平稳的,也是各态遍历的。■ 2. 性质 下面给出一些平稳随机信号的自相关函数的性质。 性质1若 X(n)是实信号,则RX(m)是偶函数,即满足 RX(m)=RX(-m)(5.2.18) 若 X(n)是复信号,则 RX(-m)=R*X(m)(5.2.19) 下面以实信号为例证明。 证明利用自相关函数的定义式(5.2.17),并令n+m=k,则n=k-m,代入公式得 RX(m)=E{X(n)X(n+m)}=E{X(n+m)X(n)}=RX(-m) 性质2RX(0)≥RX(m) 即当 m=0时,平稳过程的相关函数具有最大值,其物理意义是,同一时刻随机过程自身的相关性最强。 性质3周期平稳过程的自相关函数必是周期函数,且与过程的周期相同。 若平稳过程 X(n)满足条件 X(n)=X(n+N),则称它为周期平稳过程,其中N为过程的周期。 性质4RX(0)=E{X2(n)} 即平稳过程的均方值可以令自相关函数中 m =0得到。RX(0)代表了平稳过程的“总平均功率”。 性质5不包含任何周期分量的非周期平稳过程满足 limm→∞RX(m)=RX(∞)=μ2X(5.2.20) 从物理意义上讲,当 m 增大时,X(n)与 X(n+m)之间的相关性会减弱,在 m→∞ 的极限情况下,两者互相独立,于是有 limm→∞RX(m)= limm→∞E{X(n)X(n+m)}= limm→∞E{X(n)}E{X(n+m)}=μ2X 故有 RX(∞)=μ2X 3. 估计方法 在实际应用中,我们所遇到的信号大都是实际的物理信号,因此是因果性的,即当 n<0 时,x(n)≡0,且x(n)是实信号,其自相关为 RX(m)= limN→∞1N∑N-1n=0x(n)x(n+m)(5.2.21) 在实验中,只能得到N个观测值。如何由这N个值估计出x(n)的自相关函数?下面介绍自相关函数估计的直接法和快速计算法。 1) 自相关函数的直接估计 若观测数据的点数N为有限值,根据式(5.2.21),估计RX(m)的一种方法是 R^x(m)=1N∑N-1n=0x(n)x(n+m) 由于x(n)只有N个观测值,因此对于每一个固定的延迟 m,可以利用的数据只有(N-1-|m|)个,在 0~N-1 的范围内,上式变为 R^x(m)=1N∑N-1-|m|n=0x(n)x(n+|m|)(5.2.22) R^x(m)的长度为2N-1,它是以 m=0 为偶对称的。对于一个固定的 N,只有当|m|N时,R^x(m)的均值才接近真值RX(m); 当|m| 越接近于N时,估计的偏差越大。 对RX(m)的另一种直接估计方法是对式(5.2.22)稍做修改,即为 R^x(m)=1N-|m|∑N-1-|m|n=0x(n)x(n+m)(5.2.23) 该式是对RX(m)的无偏估计,也是一致估计。 2) 自相关函数的快速计算 利用式(5.2.22)和式(5.2.23)计算R^x(m)时,如果 m 和N都比较大,则需要的乘法次数太多,运算量太大,因此,其应用受到限制。可以利用FFT来实现R^x(m)的快速计算,式(5.2.20)可以写成 R^x(m)=1N∑N-1n=0xN(n)xN(n+m)(5.2.24) 式中,xN(n)表示由在x(n)中的N个观测值所组成的序列。对R^x(m)求傅里叶变换得 ∑N-1m=-(N-1)R^x(m)e-jωm=1N∑N-1m=-(N-1)∑N-1n=0xN(n)xN(n+m)e-jωm =1N∑N-1n=0xN(n)∑N-1m=-(N-1)xN(n+m)e-jωm 因为两个长度为N的序列的线性卷积,其结果是一个长度为2N-1点的序列,因此,为了能用DFT来计算线性卷积,需要把这两个序列的长度扩充到2N-1点,利用DFT计算相关时,也是如此。为此,把 xN(n)补N个零,得 x2N(n),即有 x2N(n)=xN(n),n=0,1,…,N-1 0,N≤n≤2N-1 (5.2.25) 记 x2N(n)的傅里叶变换为 X2N(ejω),则 ∑N-1m=-(N-1)R^x(m)e-jωm=1N∑2N-1n=0x2N(n) ejωn∑N-1m=-(N-1)x2N(n+m)e-jω(n+m)(5.2.26) 令 l=n+m,由于 x2N(n+m)=x2N(l)的取值范围是 0~2N-1,所以,l 的变化范围也应是 0~2N-1。这样,式(5.2.26)右边为 1N∑2N-1n=0x2N(n)ejωn∑2N-1l=0x2N(l)e-jωl=1NX2N(ejω)2 所以, ∑N-1m=-(N-1)R^x(m)e-jωm=1NX2N(ejω)2(5.2.27) 式中,X2N(ejω)2 是有限长信号 x2N(n)的能量谱,除以N后即为功率谱。这说明,由式(5.2.22)估计出的自相关函数R^x(m)和 x2N(n)的功率谱是一对傅里叶变换。X2N(ejω)可以用FFT快速计算,所以,用FFT计算自相关函数的一般步骤[3]如下: ① 对 xN(n)补N个0,得 x2N(n),对 x2N(n)进行DFT,得 X2N(k),0,1,…,2N-1; ② 求 X2N(k)的幅值平方,然后,除以 N,得 1NX2N(k)2; ③ 对 1NX2N(k)2 作逆变换,得R^0(m); ④ 将R^0(m)中的N≤m≤2N-1的部分向左平移2N点后形成R^x(m)。 4. 自相关函数的应用 1) 检测淹没在随机噪声中的周期信号 设正弦信号x(t)=x0sin(Ωt+),则有 Rx(τ)= limT→∞1T∫T/2-T/2x20sin(Ωt+)sin[Ω(t+τ)+]dt(5.2.28) 令Ωt+=α,则dt=1Ωdα,且ΩT=2π,故得 Rx(τ)=x202π∫π-πsinα[sinαcos(Ωτ)+sin(Ωτ)cosα]dα=x202cos(Ωτ) (5.2.29) 可见,正弦信号的自相关函数是余弦函数,两者频率相同,而随机噪声随时间滞后增加,其前后相似程度迅速减弱,自相关函数趋于0。因此,在一定延时后,自相关函数的周期性反映了原信号中含有的周期成分。 例如,在水域中探索有无潜艇通过。潜水艇的发动机发出的是周期性信号,而海浪是随机的,如果经过相关分析发现有周期性峰值,就可以得知可能有潜艇通过。又如,汽车在砂石路上行驶时,测出车桥上的振动加速度十分杂乱。但是,通过自相关分析,可以看出车桥的振动有一定的周期性。图5.2.1是北京越野车在某砂石路上以20km/h车速行驶的试验结果。从图中可以看出,延迟0.1s就有峰值,汽车车速20km/h(5.55m/s),而在路面上大概相隔0.55m就有一个起伏。汽车行驶时,1s中大概通过了10个起伏。自相关函数的分析结果证实了这一点[4]。 图5.2.1汽车前桥和车身上振动加速度自相关函数 (采样间隔为5ms,采样10段数据,压电式加速度传感器) 2) 检测信号的回声 若信号中存在有时间延时 τ0的回声或反射,那么自相关函数将在 τ=τ0处达到峰值,而且归一化自相关函数Rx(τ)/Rx(0)的大小将给出回声的相对强度的量度。当信号是宽频带信号或在空中发射的雷达信号时,经目标反射产生回波,对回波做相关分析,就可确定目标的距离、方位和速度。 5.2.2随机信号的互相关函数及其应用 1. 定义 1) 一般随机信号 互相关函数是用来描述两个随机过程 X(t)和 Y(t)之间统计关联特性的数字特征,其定义为 RXY(t1,t2)=E{X(t1)Y(t2)}=∫∞-∞∫∞-∞xypX,Y(x,y;t1,t2)dxdy(5.2.30) 式中,pX,Y(x,y;t1,t2)是两个随机过程 X(t)和 Y(t)的联合概率密度。 若X(n)和 Y(n)为两个离散随机信号,则互相关函数定义为 RXY(n1,n2)=E{X*(n1)Y(n2)} (5.2.31) 2) 广义平稳随机信号 设两个广义平稳随机信号 X(t)和 Y(t),则其互相关函数为 RXY(t1,t2)=E{X(t1)Y(t2)}=RXY(τ),τ=t2-t1 若是离散信号 X(n)和 Y(n),则互相关函数为 RXY(n1,n2)=E{X*(n)Y(n+m)}(5.2.32) 3) 各态遍历随机信号 互相关函数为 Rxy(τ)=E[x(t)y(t+τ)]= limT→∞12T∫T-Tx(t)y(t+τ)dt(5.2.33) Rxy(m)=E{x(n)y(n+m)}= limN→∞12N+1∑Nn=-Nx(n)y(n+m)(5.2.34) 例5.2.3从含有噪声的记录中检查信号的有无。 设一个随机信号x(n)中含有加性噪声 u(n),并且可能含有某个已知其先验知识的有用信号 s(n),即 x(n)=s(n)+u(n) 为了检查x(n)中是否含有 s(n),可以对x(n)和 s(n)作互相关,因此有 Rsx(m)=E{s(n)x(n+m)}=E{s(n)s(n+m)+s(n)u(n+m)} 一般认为信号和噪声是不相关的,即 E{s(n)u(n+m)}=0,所以 Rsx(m)=E{s(n)s(n+m)}=Rs(m) 这样,可以根据作互相关的结果是否与Rs(m)相符合来判断x(n)中是否含有s(n)。 ■ 例5.2.4测定系统的频率响应。 为了测定一个未知参数的线性系统的频率响应,对它输入一个功率为1的白噪声序列 u(n),记其输出为 y(n),计算输入和输出的互相关为 Ruy(m)=E{u(n)y(n+m)}=Ru(m)*h(n) 因为Ru(m)为一δ函数,所以Ruy(m)=h(n)。 对Ruy(m)作傅里叶变换,便可以得到 Puy(ejω)=H(ejω)。 ■ 2. 性质 互相关函数具有下列性质。 (1) 互相关函数与均值μ、标准差σ有如下关系: -σxσy+μxμy≤Rxy(τ)≤σxσy+μxμy(5.2.35) 式中,μx,μy 和σx,σy 分别是x(t),y(t)的均值和标准差,如图5.2.2所示[1]。 图5.2.2两个平稳随机过程x(t)和 y(t)的互相关函数的性质图示 (2) Rxy(τ)不是偶函数,是不对称的,与自相关函数不同。另外, Rxy(τ)=Ryx(-τ) (3) 如果x(t)与 y(t)是两个完全独立无关的零均值信号,则Rxy(τ)=0,所以互相关函数能够捡拾隐藏在外界噪声中的规律性信号。 (4) Rxy(τ)的最大峰值一般不在 τ=0 处。图5.2.3是一随机时间历程记录互相关函数Rxy(τ)与时间位移 τ之间的关系。如果x(t)是对一系统的输入信号,而y(t)是系统的输出信号,则由最高峰处读出的τ就是该系统的滞后时间。互相关函数的这一主要特性在工程上得到了广泛的应用。 图5.2.3典型互相关图 (5) 与自相关分析一样,在进行互相关分析时,关键问题是选择 Δt,最好对峰值出现的位置要有估计,使之不要出现在互相关图以外,当然,也不要过分靠近纵轴线,这样测出的 τ值精度不高。一般可以先选较大的Ts做一次,以便看清Rxy(τ)的全貌,然后,再选择适当的 Ts 进行分析。 3. 估计方法 如同自相关函数一样,互相关函数也有两种估计方法,即直接方法和间接方法。 1) 直接方法 互相关函数的无偏估计定义为 R^xy(m)=1N-|m|∑N-1-|m|n=0x(n)y(n+m)(5.2.36a) R^yx(m)=1N-|m|∑N-1-|m|n=0y(n)x(n+m)(5.2.37a) 互相关函数的有偏估计定义为 R^xy(m)=1N∑N-1-|m|n=0x(n)y(n+m)(5.2.36b) R^yxm=1N∑N-1-|m|n=0ynxn+m(5.2.37b) 式(5.2.36a)~式 (5.2.37b)表明,随着m的增加,互相关分析中的相乘相加的数据量逐渐减少。但是,在有偏估计中,分母N保持不变; 而在无偏估计中,分母N-m也在同步地减小。因此,互相关分析的有偏估计和无偏估计存在明显的差异,即在互相关分析的有偏估计中,随着m的增加,理论上会对互相关分析结果的幅值进行衰减,且m越大,衰减越大。而在互相关分析的无偏估计中,并不会对互相关分析结果的幅值进行衰减。互相关分析的有偏估计和无偏估计的这一特点,使得其应用在延迟时间T的计算中具有不同的适用范围[5]。 (1) 当对周期信号进行互相关分析时,在有偏估计中,随着m的增加,理论上会对互相关分析结果的幅值进行衰减,且m越大,衰减越大。此外,周期信号的互相关分析结果为同周期的周期信号。因此,在有偏估计的结果中,存在多个相距Tm(Tm为信号周期)的局部峰值点,且局部峰值点的幅值逐渐衰减。这样,第一个局部峰值点即为互相关分析结果的峰值点。因此,要想通过找到互相关分析结果的峰值点所对应的时间,来计算两个周期信号之间的延迟时间T,就必须保证τ小于Tm,以保证τ所对应的局部峰值点为第一个局部峰值点,即互相关分析结果的峰值点。而在无偏估计中,由于理论上并不会对互相关分析结果的幅值进行衰减。因此,周期信号的互相关分析结果仍为同周期的周期信号,其各局部峰值点的幅值相等,没有单个的峰值点,从而无法通过找互相关分析结果中峰值点的方法来计算延迟时间τ。 (2) 当对局部峰值点和周期变化较小的一类准周期信号进行互相关分析时,由于这一类准周期信号与周期信号非常相似,只是局部峰值点和周期存在较小的变化。因此,与采用有偏估计分析周期信号的结果一样,存在多个局部峰值点,且对局部峰值点的幅值进行衰减,m越大,衰减越大。为了能够通过寻找互相关分析结果的峰值点来计算延迟时间τ,要求延迟时间τ小于信号的近似周期Tm,以保证第一个局部峰值点所对应的时间即为延迟时间τ。 对于无偏估计,由于这一类准周期信号的周期不完全相等,各局部峰值点也不完全相同。因此,其互相关分析结果中的各局部峰值点的幅值不完全相等。所以,无论延迟时间τ是否小于信号的近似周期Tm,只要当m的取值所对应的时间等于延迟时间τ时,互相关分析结果的幅值最大。所以,可以通过寻找互相关分析结果峰值点来计算延迟时间τ。但是,该方法要求两信号完全相似,即两个信号之间除了存在一定的相位差外,只有幅值上的按比例衰减。这是因为无偏估计不对互相关分析结果进行衰减,导致其结果中虽然峰值点所对应的时间为延迟时间τ,但是,各局部峰值点的幅值基本相同。若有干扰存在,就将导致其中一路信号受到干扰,可能造成互相关分析结果中某个局部峰值点的幅值变大,超过延迟时间τ所对应的峰值点的幅值。此时,通过寻找互相关分析结果的峰值点来计算延迟时间τ,就会产生较大的误差。由于在实际应用中很难保证两路输出信号完全相似,所以,这种方法的抗干扰能力较差。 (3) 当对冲激信号进行互相关分析时,由于冲激信号在大部分时间内幅值为0(或白噪声),其有用信号仅占其中一小部分时,只有当m的取值所对应的时间等于两信号间的延迟时间τ时,非零点的相乘相加量最多,互相关分析结果取最大值。因此,无论采用有偏估计还是无偏估计,均可以通过寻找互相关分析结果峰值点来计算延迟时间τ,且对延迟时间τ的大小没有要求。 2) 间接方法 假定x(n)与 y(n)的初始采样容量为 N=2p,用FFT方法计算互相关函数时的过程是: 先通过FFT求得互功率谱函数,然后计算互功率谱函数的逆傅里叶变换,它包括两个独立的FFT运算,一个是对x(n)的计算,另一个是对 y(n)的计算; 然后计算FFT逆变换,就可得到互相关函数。 4. 互相关函数的应用 互相关分析的主要应用如下。 1) 测量滞后时间 若系统是线性的,则计算输入和输出之间的互相关,就可以得到滞后时间。当系统的输出与输入之间有一定的时间差时,互相关函数在时间差等于信号通过系统所需时间值时将出现峰值。实际上,互为线性关系的两个信号,其平均乘积在信号间的时间差为零时总是取得最大值。因此, 可以采用输入输出互相关图中峰值的时间差来确定系统的时间滞后。 当信号由a点传播到b点时,两点上的信号 xa(t)与 xb(t)之间的互相关函数Rab(τ)将在 τ=τL 处出现峰值。 设信号传播速度为 V,a和b两点距离 L,则信号由a点传播到b点的时间延迟 τL 为 τL=LV (5.2.38) 式(5.2.38)中的三个参数,已知任两个参数,便可确定第三个参数,这在工程实践中有很多应用。图5.2.4是测定热轧钢运动速度的示意图。已知两光电接收器的距离为d,则可得出两个光电信号的互相关峰值处的时延 τd,因此,钢带运动速度为 V=d/τd 图5.2.4热轧钢运动速度的测定示意图 又如,微血管中红细胞运动速度的测量,见图5.2.5[6]。利用微循环血流仪可以方便地把手指尖部微血管中红细胞运动的状态显示在屏幕上。为了测量红细胞的运动速度,在屏幕上所显示的微血管的上、下游设置了两个窗口 W1 和 W2,当红细胞经过这两个窗口时,将产生两个光密度信号 x1(t)和 x2(t)。若假定红细胞在血管中 图5.2.5微血管中红细胞运动速度的测量 运动时不变形,在其他情况下亦保持不变,则 x1(t)和 x2(t)应完全相似。由于信号是在一段有限长时间内提取的,可以计算其互相关函数为 R12(τ)=∫∞-∞x1(t)x2(t+τ)dt(5.2.39) R12(τ)取得极大值时的延迟 τ0便是红细胞从窗口 W1 运动到 W2 时所需要的时间,求出 τ0,便可计算出红细胞的运动速度 v,即有 v=L/τ0(5.2.40) 式中,L是窗口 W1 和 W2 之间的距离,可以事先测出。 研究一个控制系统,其系统的滞后时间是一个很重要的参数。例如,进行汽车操纵性试验时,汽车操纵系统的反应时间是一个重要的测定参数,但采用的方法由于起点不好确定及人的主观因素,常常可能有1/3的误差。一般做这种方向盘阶跃响应试验,危险性较大,并且需要很大的场地。 图5.2.6汽车脉冲试验互相关函数 (车速为60km/h,采样间隔为30ms) 改用互相关分析方法来测定时,可在比较宽的路段上,固定车速直线行驶,司机间隔地、脉冲地转动方向盘使汽车车身在行驶方向上产生晃动。记录下司机转动方向盘的转角作为系统的输入,汽车垂直轴晃动的角速度作为系统的输出,进行互相关分析,得出的结果如图5.2.6所示。峰值偏离轴线的距离τ,即为汽车操纵系统的滞后时间,可以用于评价汽车执行司机转动方向盘指令的反应快慢。用互相关分析方法,可以排除人为因素的影响,精度较高,实验安全,不需要很大的场地。 用互相关方法测定深埋地下的输油管和水管的裂损位置,见图5.2.7。由于漏损处液体流动发生声源向两端传播,因此,在管道两端上分别放有传感器 x 和 y。因为漏损处离两端点距离不等,因此,声响传到传感器处有时差。在互相关图上分析出其时差为 τ′,则漏损处离两端点的中心线距离 s=12vτ′,其中,v 为声响通过管道的传播速度。如果 v 为未知,也可在一测点处敲击,用互相关分析就可以测出声响从一端到另一端所需的时间,如果测得 τ′=0,则说明漏损处正好是两测点的中点。用这方法可以避免全线开挖,从而节省了巨大的工作量。 图5.2.7用互相关法测定油管裂损位置示意图 2) 确定传递通道 互相关函数可用来确定传递通道。例如,重型机械的运转,常会引起噪声和振动,其能量可由几个通道建筑物结构传递,在有效地控制噪声和振动之前,必须精确地测出传递通道,这就可用系统输入和输出之间的互相关性来解决。由于每一条通道一般都具有不同的滞后时间,因此,在互相关图中将出现各自的峰值。如果通过另外的方法计算出各个通道的滞后时间,则这些滞后时间就可以与从互相关中峰值得到的滞后时间进行对比,从而就可确定对输出有显著影响的通道。 又如,测试车辆振动传递通道,其测试框图如图5.2.8所示,用以检查汽车司机座位的振动是由发动机引起的,还是由车轮引起的。测试方法是: 在发动机、司机座位和后轮轴上布置加速度计,经分析,发现发动机与司机座位之间的相关性较弱,而司机座位与后轮之间的相关函数则出现明显的相关。因此可以认为,司机座位的振动主要由后轮传递来的[1]。 图5.2.8车辆振动传递通道的测试框图 3) 检测噪声中的信号 互相关函数也可以用来检测隐藏在外界噪声中的信号,而且信号可以不是周期形式的。由于自相关分析不能从外界噪声中分离出随机信号,因此,需要用互相关函数。具体地讲,假定已有一个需要测定的无噪声信号(随机的或周期的),则从信号和噪声的互相关性可得到信号的互相关函数。另外,对于周期信号,在任一给定的输入信噪比和样本记录长度的条件下,由互相关函数得到的输出信噪比,要比用自相关函数得到的输出的信噪比大。 4) 系统识别 所谓系统识别就是根据实验数据,利用相关技术计算出反映系统特性的频率响应函数,即由互相关函数的傅里叶变换,求得系统的频率响应或互谱密度函数。 为了测定一个未知参数的线性系统的频率响应,对它输入一白噪声序列x(n),记录其输出 y(n),计算输入和输出的互相关为 Rxy(m)=E{x(n)y(n+m)}=Ex(n)∑∞k=-∞h(k)x(n+m-k) =∑∞k=-∞h(k)E{x(n)x(n+m-k)} =∑∞k=-∞h(k)Rxx(m-k)=Rxx(m)*h(m) (5.2.41) 因为x(n)的自相关函数Rxx(m)为一δ函数,所以 Rxy(m)=h(m)(5.2.42) 对 h(m)作离散傅里叶变换,便可求出系统的频率响应 H(k)。该过程的原理框图如图5.2.9所示。 图5.2.9测定系统的频率响应原理框图 5.3随机信号的功率谱估计 对于确定性信号,如果在时域进行分析比较复杂,则可以利用傅里叶变换将其转换到频域进行分析。同样,对于随机过程,也可以利用傅里叶变换来分析其频谱结构,不过,随机过程的样本函数一般不满足傅里叶变换的绝对可积条件 ∫+∞-∞x(t)dt<∞ (5.3.1) 而且随机过程的样本函数往往不具有确定的形状。一般,平稳随机过程的单个样本的平均功率是有限的,即 Gξ= limT→∞12T∫T-Tx(t)2dt<∞ (5.3.2) 式中,ξ 表示某个实验结果。若x(t)为随机过程 X(t)的样本函数,X(t)代表噪声电流或电压,则 Gξ 表示x(t)消耗在1Ω电阻上的平均功率。这样,对于随机过程的样本函数而言,研究它的频谱没有意义,研究其平均功率随频率的分布才有意义。 5.3.1随机信号的功率谱密度 首先把随机过程 X(t)的样本函数x(t)任意截取一段,长度为2T,记作 xT(t), xT(t)=x(t),t≤T 0,t>T(5.3.3) 称为截断函数,如图5.3.1所示。 图5.3.1x(t)及其截断函数 对于有限持续时间的 xT(t)而言,傅里叶变换是存在的,所以有 XT(Ω)=∫∞-∞xT(t)e-jΩtdt =∫T-TxT(t)e-jΩtdt (5.3.4) xT(t)=12π∫∞-∞XT(Ω)ejΩtdΩ(5.3.5) XT(Ω)即为 xT(t)的频谱函数[2]。 由于x(t)是随机过程 X(t)的一个样本函数,取哪一个样本函数取决于实验结果 ξ,且是随机的,因此XT(Ω)和 xT(t)也都是实验结果 ξ 的随机函数,最好写成 XT(Ω,ξ)和 xT(t,ξ)。 现将式(5.3.5)代入式(5.3.2),并考虑到讨论的是实随机过程,x(t)是实函数,可得某个样本函数的平均功率 Gξ 为 Gξ= limT→∞12T∫T-TxT(t,ξ)2dt = limT→∞12T∫T-TxT(t,ξ)12π∫∞-∞XT(Ω,ξ)ejΩtdΩdt = limT→∞12T∫T-T12πXT(Ω,ξ)∫T-TxT(t,ξ)ejΩtdtdΩ = limT→∞12T∫∞-∞12πXT(Ω,ξ)2dΩ =12π∫∞-∞limT→∞12TXT(Ω,ξ)2dΩ(5.3.6) 如果频率函数满足以下两个条件,则为信号的功率谱密度函数。一是当在整个频率范围内对它进行积分以后,就给出信号的总功率; 二是它描述了信号功率在各个不同频率上分布的情况。式(5.3.6)的被积函数 limT→∞12TXT(Ω,ξ)2 具备了上述特性,它表示随机过程的某一个样本函数 x(t,ξ)在单位频带内、消耗在1Ω电阻上的平均功率,因此,称它为样本函数的功率谱密度函数,记作 PX(Ω,ξ)。 PX(Ω,ξ)= limT→∞12TXT(Ω,ξ)2(5.3.7) 如果对所有实验结果 ξ 取统计平均,得 PX(Ω)=EPX(Ω,ξ)=ElimT→∞12TXT(Ω,ξ)2= limT→∞12TEXT(Ω,ξ)2 (5.3.8) 这里的 PX(Ω)是Ω的确定函数,不再具有随机性,它表示随机过程 X(t)在单位频带内在1Ω电阻上消耗的平均功率。因此,PX(Ω)被称为随机过程 X(t)的功率谱密度函数,简称功率谱密度。 如果将式(5.3.6)对所有的实验结果 ξ 取统计平均,则得到随机过程 X(t)的平均功率为 G=E{Gξ}= limT→∞12T∫T-TEx(t,ξ)2dt= limT→∞12T∫T-TEX(t)2dt =12π∫∞-∞limT→∞12TEXT(Ω,ξ)2dΩ=12π∫∞-∞PX(Ω)dΩ(5.3.9) 可见,随机过程的平均功率可以由它的均方值的时间平均得到,也可以由它的功率谱密度在整个频率域上积分得到。 若 X(t)为平稳过程时,此时均方值为常数,于是,式(5.3.9)可以写成 G=E{Gξ}=E{X2(t)}=12π∫∞-∞limT→∞12TEXT(Ω)2dΩ(5.3.10) 或 G=E{X2(t)}=RX(0)=12π∫∞-∞PX(Ω)dΩ(5.3.11) 该式说明,平稳过程的平均功率等于该过程的均方值,它可以由随机过程的功率谱密度在全频域上的积分得到。 若 X(t)为各态历经过程,则有 PX(Ω)= limT→∞12TXT(Ω,ξ)2 (5.3.12) PX(Ω)是从频率角度描述 X(t)的统计特性的重要数字特征的,但是,PX(Ω)仅表示 X(t)的平均功率在频域上的分布情况,不包含 X(t)的相位信息。 例5.3.1随机过程 X(t)为 X(t)=acos(Ω0t+Φ) 式中,a 和Ω0是常数; Φ是在(0,π/2 )上均匀分布的随机变量,求随机过程 X(t)的平均功率G。 解 E{X2(t)}=E{a2cos2(Ω0t+Φ)}=Ea22+a22cos(2Ω0t+2Φ) =a22+a22∫π/202πcos(2Ω0t+2φ)dφ=a22-a2πsin(2Ω0t) 显然,这个过程不是平稳过程,所以,必须做一次时间平均,由下式求其平均功率,得 G= limT→∞12T∫T-TE{X(t)2}dt= limT→∞12T∫T-Ta22-a2πsin(2Ω0t)dt=a22 5.3.2功率谱密度的性质 功率谱密度有如下的性质。 性质1非负性,即 PX(Ω)≥0 (5.3.13) 根据功率谱密度的定义式(5.3.8),考虑到其中的 XT(Ω)2 必然非负,故其数学期望值也是非负的,因而得证。 性质2 PX(Ω)是实函数。 因为 XT(Ω)2 是实函数,所以它的数学期望必然是实的。 性质3当随机信号是实过程时,其功率谱 PX(Ω)是偶函数,即 PX(Ω)=PX(-Ω)(5.3.14) 因为 XT(Ω)2=XT(Ω)XT(-Ω),所以 XT(Ω)2 是偶函数,故它的数学期望也必然是偶函数。 5.3.3功率谱密度与自相关函数的关系 功率谱密度从频率角度描述过程统计特性的数字特征,而相关函数则从时间角度描述过程统计特性的最主要的数字特征,两者描述的是同一个对象,因此,它们之间必然有一定的关系。下面将证明,平稳过程在一定的条件下,其自相关函数RX(τ)和功率谱密度 PX(Ω)构成傅里叶变换对。 证明从式(5.3.8)出发,并考虑到 XT(Ω,ξ)=∫∞-∞xT(t,ξ)e-jΩtdt=∫T-TxT(t,ξ)e-jΩtdt XT(Ω,ξ)2=XT(Ω,ξ)XT(-Ω,ξ) (5.3.15) 于是,式(5.3.8)可以写成 PX(Ω)= limT→∞E12T ∫T-TxT(t1,ξ)ejΩt1dt1∫T-TxT(t2,ξ)e-jΩt2dt2 =limT→∞12T ∫T-T∫T-TE[XT(t1)XT(t2)]e-jΩ(t2-t1)dt1dt2 (5.3.16) 现交换积分次序,并引入下面相关函数的表示: RXT(t1,t2)=EXT(t1)XT(t2),-T<(t1,t2)<T (5.3.17) 注意,有上述时间限制,式(5.3.17)与 X(t)的相关函数是一样的,则式(5.3.16)改写为 PX(Ω)= limT→∞12T∫T-T∫T-TRX(t1,t2)e-jΩ(t2-t1)dt1dt2 (5.3.18) 在式中作积分变量替换,则有 t=t1,dt=dt1 τ=t2-t1=t2-t,dτ=dt2 于是,式(5.3.18)可改写为 PX(Ω)= limT→∞12T∫ T-t-T-t∫T-TRX(t,t+τ)dte-jΩτdτ 将极限符号写入,得 PX(Ω)=∫∞-∞limT→∞12T∫T-TRX(t,t+τ)dte-jΩτdτ(5.3.19) 其中,括号里的量可以看作是非平稳过程自相关函数的时间平均,即为 RX(τ)= limT→∞12T∫T-TRX(t,t+τ)dt (5.3.20) 于是,式(5.3.19)改写为 PX(Ω)=∫∞-∞ RX(τ)e-jΩτdτ (5.3.21) 即时间平均自相关函数与功率谱密度为傅里叶变换对。现假设 X(t)为平稳过程,则时间平均自相关函数等于集合平均自相关函数,式(5.3.19)为 PX(Ω)=∫∞-∞RX(τ)e-jΩτdτ(5.3.22) 可见,平稳过程的功率谱密度就是其自相关函数的傅里叶变换。若进行傅里叶反变换,则有 RX(τ)=12π∫∞-∞P(Ω)ejΩτdΩ(5.3.23) 证毕。 式(5.3.22)和式(5.3.23)是平稳过程统计特性频域描述(功率谱密度)和时域描述(相关函数)之间的重要关系式。这对关系式在实际中有着广泛的应用价值。由于这对关系式是由美国学者维纳(Wiener)和苏联学者辛钦(χИНЧин)得出的,因此也常被称为维纳辛钦定理或维纳辛钦公式。 在式(5.3.23)中,当 τ=0,则可得式(5.3.11)。实际上,当满足∫∞-∞P(Ω)dΩ<∞,或者说平均功率有限时,式(5.3.23)才能成立。这个条件在实际应用中通常是满足的。 在以上的介绍中,PX(Ω)应分布在-∞~+∞ 的频率范围内,实际上,负频率是不存在的。在公式中,频率从负到正,纯粹只有数学上的意义和为了运算方便。有时也采用另一种功率谱密度,即“单边”谱密度,也称作“物理”功率谱密度,记作 FX(Ω)。FX(Ω)只分布在Ω≥0 的频率范围内,FX(Ω)与 PX(Ω)的关系是 FX(Ω)=2PX(Ω),Ω≥0 0,Ω<0(5.3.24) 还需要指出的是,以上所介绍的功率谱密度都属于连续的情况,这意味着相应的随机过程不能含有直流成分或者周期性成分,这也是式(5.3.22)傅里叶变换要求RX(τ)绝对可积的条件。这是因为,功率谱密度是指“单位带宽上的平均功率”,而任何直流分量和周期性分量,在频域上都表现为频率轴上某点的零带宽内的有限平均功率,都会在频域的相应位置上产生离散谱线,而且在零带宽上的有限功率等效于无限的功率谱密度。于是,当平稳包含有直流成分时,其功率谱密度在零频率上应是无限的,而在其他频率上是有限的。换句话说,该过程的功率谱密度函数曲线将在相应的零频率点上存在δ函数。同理,若平稳过程含有某个周期成分,则其功率谱密度函数曲线将在相应的离散频率点上存在δ函数。借助于δ函数,维纳辛钦公式就可推广应用到这种含有直流或周期性成分的平稳过程中来。具体来说,一是当平稳过程有非零均值时,正常意义下的傅里叶变换不存在。但是,非零均值可以用频域原点处的δ函数表示,该δ函数的权重即为直流分量的功率。二是当平稳过程含有对应于离散频率的周期分量时,该成分就在频域的相应频率上产生δ函数。 δ函数的基本而重要的性质是筛选特性,即对任一连续函数 f(t)有 ∫∞0δ(t)f(t-τ)dt=f(τ)(5.3.25) 由此可以写出下面重要的傅里叶变换对 12π∫∞-∞e-jΩτdτ=δ(Ω)12π=12π∫∞-∞δ(Ω)ejΩτdΩ(5.3.26) ∫∞-∞δ(τ)ejΩτdτ=1δ(τ)=12π∫∞-∞e-jΩtdΩ(5.3.27) 例5.3.2已知随机过程 X(t)的自相关函数为 RX(τ)=12cos(Ω0τ) 求功率谱密度。 解 PX(Ω)=∫∞-∞12e-jΩτcos(Ω0τ)dτ=12∫∞-∞12(ejΩ0τ+e-jΩ0τ)e-jΩτdτ =14∫∞-∞ej(Ω-Ω0)τdτ+14∫∞-∞ej(Ω+Ω0)τdτ=π2δ(Ω-Ω0)+δ(Ω+Ω0) PX(Ω)的图形如图5.3.2所示。 ■ 图5.3.2例5.3.2的 PX(Ω)的图形 上述维纳辛钦定理不仅适用于连续时间随机信号,而且也适用于离散时间随机序列,即对具有零均值实平稳随机序列的功率谱密度 PX(ω)与序列的自相关函数RX(m)是一对离散时间傅里叶变换对。 设 X(n)为一具有零均值的平稳随机序列信号,其自相关函数为 RX(m)=EX(n)X(n+m) (5.3.28) 当RX(m)满足绝对可积条件时,把均匀间隔的离散时间信号的自相关函数作离散时间傅里叶变换,便得到其功率谱为 PX(ω)=∑∞m=-∞RX(m)e-jmω(5.3.29) 同时,有 RX(m)=12π∫π-πPX(ω)ejmωdω (5.3.30) 离散时间信号的功率谱主要有以下特点。 (1) 功率谱是周期性的,因此可以作傅里叶级数分解,而RX(m)正是分解出的各次谐波的系数。 (2) 反演变换的积分区间是[-π,π]。 计算式(5.3.29)时,可以先对RX(m)两边做z变换,得 PX(z)=∑∞m=-∞RX(m)z-m(5.3.31) 再令 z=ejω,便得到功率谱。 例5.3.3设RX(m)=am,a<1,m=0,±1,±2,…,求功率谱 PX(ω)。 解根据式(5.3.31),有 PX(z)=∑∞m=-∞RX(m)z-m=∑∞m=0amz-m+∑-∞m=-1a-mz-m=11-az-1+∑∞m′=1am′zm′ =11-az-1+az1-az=1-a21+a2-a(z-1+z) 因此,其功率谱为 PX(ω)=1-a21+a2-2acosω ■ 注意,为了书写的方便,我们有时把 ejω 的函数写成ω 的函数,例如将 PX(ejω)写成 PX(ω),将P^X(ejω)写成 P^X(ω)。 5.3.4功率谱估计的方法 在生产和工程中,人们所要分析和处理的实际信号往往具有随机性,即使是确定信号,在传输的过程中也不可避免地会受到噪声的干扰。为了排除噪声的干扰,提取有用的信号,利用在有限的时间范围内的观测数据来计算离散随机信号的统计特性。但是,由于观测数据的长度N只能取有限值,因此,所计算出的数字特征只能是 N→∞时所得的相应真值的估计。功率谱的估计在随机信号处理中有着极其广泛的应用,其主要目的是用来揭示一些看来杂乱无章、无规律可循的事物中所蕴含的周期性。目前,功率谱估计的方法很多,大致可分为以下两种。 (1) 经典法(线性估计法): 即用传统的傅里叶变换分析方法求功率谱的方法。它又可分为间接法(相关估计法)和直接法(周期图法)。间接法是由数据的自相关序列求功率谱的方法; 直接法是由数据直接用离散傅里叶变换求功率谱的方法。 (2) 现代法(非线性估计法): 即用参量信号模型来估计功率谱的方法。参量模型又可分为自回归信号模型(AR模型)、滑动平均模型(MA模型)和自回归滑动平均模型(ARMA模型)。 本书仅介绍两种经典功率谱估计方法。我们知道,傅里叶变换是揭示周期性的有力工具,功率谱估计的经典法就是建立在傅里叶变换的基础上的。另外,在下面的讨论中,假设随机信号是各态历经的,所以,集合平均可以由单一样本的时间平均来实现,故将原来功率谱密度和自相关函数符号中大写字母的下标换成小写字母的下标。 1. 自相关法 自相关法是根据维纳辛钦定理来计算功率谱的。这种方法是Blackman(布莱克曼)与Tukey(图基)两人于1958年提出的,所以,也称为BT法。具体做法是,先根据样本 x(m)的N个观测值 xN(m)估计出自相关函数R^x(m),然后求R^x(m)的傅里叶变换,即得 xN(m)的功率谱。现在考虑零均值平稳遍历随机信号,并利用FFT算法来计算功率谱。随机信号 xN(m)的自相关函数为 R^x(m)=1N[x(m)*x(-m)],-(N-1)≤m≤N-1 (5.3.32) 序列 x(m)的长度为N,故R^x(m)的长度为2N-1,因此,应将此线性卷积等同为长度为2N-1的循环卷积,用2N-1点的FFT算法来进行计算功率谱。 R^x(m)算出后,根据离散信号的傅里叶变换式即可求得功率谱估计为 P^(ω)=∑∞m=-∞R^x(m)e-jωm=∑N-1m=-(N-1)R^x(m)e-jωm (5.3.33) 将数字频率ω 离散化得 P^x(k)= P^x(ω)|ω=2πNk=∑N-1m=-(N-1)R^x(m)e-j2πNkm(5.3.34) 2. 周期图法 周期图法是把随机信号x(n)的N点观测数据 xN(n)视为一能量有限信号,对其直接进行傅里叶变换,得到 XN(ejω),然后对谱的模进行平方运算求得功率谱,即 P^x(ejω)=1NXN(ejω)2 (5.3.35) 将ω 在单位圆上等间隔取值,得 P^x(k)=1NXN(k)2,k=0,1,…,N-1 (5.3.36) 由于 XN(k)可以用FFT快速计算,所以,P^x(k)也可被方便地求出。周期图这一概念是由Schuster于1899年首先提出的,因为它是直接由傅里叶变换得到的,所以,人们习惯上称为直接法。它是一种非常古老的方法,被广泛采用。在1958年维纳辛钦定理被提出之后,BT法由于算法简捷而取代了周期图法。但是,在1965年CoolyTukey(库利图基)创造性地提出了FFT算法后,周期图法又成为应用最为广泛的功率谱估计方法。 周期图法虽然简便,但是,它只是随机序列功率谱密度的估计值,因此,与功率谱密度 Px(ejω)真值之间还存在着误差。若当序列长度N趋于无限时,周期图与统计平均值的偏差或 P^x(ejω)的方差趋近于零,则认为是满意的估计。因此,实际中在信号处理技术上往往采取措施,对周期图法进行修改,以达到尽量减少估计误差(数学期望与方差)的目的,其方法有以下两种。 1) 采取窗处理减少功率泄漏 用有限长序列来估计功率谱,实质上等于加矩形窗截断随机序列,因此,必然会出现吉布斯(Gibbs)效应,使信号原来集中在小范围的功率扩散到较大的频带内。这种因截断而扩散到主瓣以外的功率称为泄漏功率,它使谱估计值与真值之间的误差加大。为此,采取窗处理的办法,先将序列x(n)乘以适当的窗函数w(n),即 XN(k)=DFT[x(n)w(n)](5.3.37) 然后,再利用FFT算法来求出相应的功率谱值。在时域加窗可以促使频域收敛得快一些,其物理含义是对功率谱给予平滑滤波,减少功率泄漏。 2) 采取平均化处理减小统计变异性 在实际中,对于平稳随机信号都是通过一次观测求出其全部统计特性。但是,这只有当观测时间趋于无限长时,才能得到正确的结果。因此,对有限长数据的处理,其统计特性一定存在着误差,这就是所谓的统计变异性。减小这种变异性的方法有两种: 一是加长观测时间,处理大量数据,其结果必然对设备提出过高的要求; 二是对同一现象独立观测多次,然后对结果进行平均化处理。例如,观察总长度为 Tl,则可把它分为等长的 K 段,每段长度为 Tr,然后利用FFT算法分别求出每段的功率谱密度估值 Pm(k),其中,m=1,2,…,K。于是,经平均化处理后的功率谱估计值为 P^x(k)=1K∑Km=1P^m(k)(5.3.38) 故得极限的频谱分辨率为 Δf=1Tr=KTl (5.3.39) 即分辨率减弱为原值的1/K。由概率论可知,K 个相同的独立随机变量的平均方差应等于单独方差的 1/K。所以,若假设各段周期图彼此互相独立无关,则经平均化处理后,将使估计的方差大约减小 K 倍。换句话说,这种方法是用降低分辨率来换取估计方差的减小的。 为了减小方差,同时考虑到分辨率的改善,Welch根据估值理论于1967年提出了一种改进的周期图平均法。该法的基本思想是,为使分段数增加而每段的长度又不减少,所以,在分段处理过程中,采用2∶1覆盖的分段。在这种方法中,由于分段的数量增加,因而方差进一步减小。例如,有一长度 N=500的数据,原分为5段,每段100个数据,现若覆盖50%,则段数将增加到10段,而每段仍为100个数据,这样既保持了一定的分辨率,又有利于估值偏差的减小。虽然由于覆盖使数据的依赖性有所增加,但是,总的方差仍然有显著的降低[8]。 可以将加窗和平均化处理技术结合起来,其求解步骤如下: (1) 将原始实序列长度为N的原始序列按2∶1覆盖分段; (2) 选用适当的时窗函数对每段进行窗处理; (3) 利用FFT算法估计每段功率谱; (4) 通过求每段功率谱之和进行平均化处理。 这种改进的周期图法的流程图如图5.3.3所示。 图5.3.3改进的周期图法流程 3. 经典谱估计法的优缺点 (1) 经典谱估计,无论是周期图法还是自相关法,都可以用FFT算法计算功率谱,且计算速度快,物理概念明确,因而仍是目前较常用的谱估计方法。 (2) 在实际应用中,存在两个不符合实际的假设: 其一,自相关法利用有限的观测数据N做自相关估计,假设N个数据之外的x(n)为0,即假设当|m|≥N时,RX(m)=0; 其二,周期图法假设数据是以N为周期的周期性延拓。这些不合理的假设,将一些不真实的信息附加在随机过程之上,限制了频率分辨率的提高和功率谱估计的质量指标,使其谱估计的均值不等于真值,而是真值功率谱与窗频谱的卷积,因而估计是有偏的; 且谱估计方差当 N→∞ 并不为零,不是一致估计。 (3) 谱的分辨率正比于 2π/N,N是所使用的数据的长度,或者说正比于被分析数据的总的时间采样长度。 (4) 由于不可避免窗函数的影响,因而使得真正的谱 P(ω)在窗口主瓣内的功率向边瓣部分“泄漏”,从而降低了分辨率。较大的边瓣有可能掩盖 P(ω)中较弱的成分,或者产生假的峰值。当分析较短数据时,这些影响更为突出。 (5) 方差性能不好,不是 P(ω)的一致估计,且N增大时,谱曲线起伏加剧。 (6) 可以通过采用不同的窗函数和频谱校正方法来改善频谱估计的性能。 5.3.5功率谱估计的应用 1. 从含有噪声的信号中确定主频率[7] 涡街流量计是一种应用比较广泛的测量气体、液体和蒸汽流量的仪表,其工作原理是卡门涡街原理。在流动的流体中放置一个阻挡物体(旋涡发生体),就会在阻挡物体的下游两侧产生两列有规律的旋涡,在满足一定条件时,产生的旋涡是稳定的,称为卡门涡街。涡街的脱离频率与流量存在对应关系。当在旋涡发生体右(或左)下方产生一个旋涡以后,在旋涡发生体上产生一个升力,在旋涡发生体内部安装的压电式传感器,将作用在旋涡发生体上的升力转换为电荷信号,电荷的变化频率与旋涡的脱离频率一致。通过检测压电传感器输出信号的变化频率,就可以得到旋涡的脱离频率,从而得到流体的流量。在理想状态下,压电式传感器的输出信号是一个正弦波,但是,在实际中,输出信号含有各种噪声,噪声包括机械振动噪声和流体流动噪声等。如何从叠加了噪声的信号中提取出涡街信号,从而提高测量频率的精度,是流量计信号处理的核心问题。现在市面上的涡街流量计均用放大、滤波、整形、计数的方法对传感器的输出信号进行处理,或是在此基础上进行数字滤波,这种基于模拟和数字电路的处理方法有某些优点,但是,这些方法从含有噪声的信号中提取信息的能力不够强,而且仪表的现场测量精度不理想。 采用周期图法直接对采样的传感器数据(也可经过加权)作傅里叶变换,再取其幅度的平方便得到功率谱。因为功率谱分析方法可将各种干扰噪声和涡街信号分离出来。因此,在信号的幅度大于干扰幅度的情况下,可以根据其能量大小判别出涡街信号的频率。 采用周期图法来进行信号处理,存在计算误差,它主要是由于数据加窗而造成的泄漏误差。在众多窗函数中,矩形窗的主瓣最窄,为了保证测量频率的精度,选用矩形窗进行截断。在周期图法进行频谱分析时,计算误差不大于频率分辨率。为此,需采用合适的频率分辨率来控制非整周期采样造成的泄漏误差。 为了提高测量精度,我们希望降低采样频率,增加采样点数。但是,采样频率的减小是有限度的,因为采样要满足采样定理。点数增加则会增大计算量,并提高对计算机数据存储量的要求。为了减少采样点数而同时满足计算精度的要求,需要分段设置采样频率。涡街信号的范围大约为2~2500Hz,采样点数为4096点,计算精度应优于0.2%,分10个频率段(或少些)设置采样频率,最低采样频率为150Hz,最高采样频率为12.4kHz。在信号频率高于30Hz时,可以满足精度要求; 当信号频率低于30Hz,如果用降低采样频率的方式来保证计算精度,采样时间会太长,此时,运用频谱分析校正方法可使计算精度优于0.2%。 2. 分析电动机噪声产生的原因[1] 电动机噪声测试框图如图5.3.4所示。被测感应电动机的极数为14级,电源频率为60Hz,电动机的转速为512r/min。噪声计用于测量电动机噪声,并与信号处理机1通道相连; 用压电式加速度计测量电动机外壳的振动加速度,经放大送入信号处理机2通道,分别计算出它们的功率谱,如图5.3.5所示,图5.3.5(a)为电动机结构加速度功率谱,图5.3.5(b)为噪声的功率谱。 图5.3.4电动机噪声测试框图 图5.3.5实测功率谱 由功率谱特性并结合电动机的转速和极数等因素综合分析,可以认为,噪声功率谱中频率为120Hz的成分是电磁噪声(为电源频率60Hz的两倍); 频率为490Hz的成分与转速和滚珠轴承的滚珠数和直径有关; 频率为1370Hz的成分也是电磁噪声。为了降低该电机的噪声声压,必须减小120Hz,490Hz和1370Hz的峰值。 3. 不解体的故障判断[4] 有许多大型设备,为了防止出现故障,需要定期检修。定期检修往往是将设备拆开来看有没有问题,这样经常拆修,机器容易损坏,而且也不能完全避免事故的发生。利用信号处理技术,可将测得的振动信号和噪声信号,用数字信号处理系统处理后,得到频谱图,根据频谱图来判断设备有无故障,以及是否需要检修。因为这种方法不必拆开设备,所以,称为不解体的故障判断,近年来得到广泛的应用,不仅应用于大型设备,如飞机、核反应堆等,也应用于发动机和齿轮箱等。图5.3.6是汽车发动机振动的频谱分析,从图中可以看到,排气门间隙过大,高频成分明显增加,因此根据高频成分大小就可以判断气门间隙调整得是否合适。 图5.3.6发动机振动的频谱图 图5.3.7是实验测得的汽车变速器的振动加速度和噪声,经过信号处理后得到的频谱图。从图中可以对比出,不正常的变速器在9.2Hz和18.4Hz 上增加了峰值,由此可以判断出变速器的某对齿轮出了故障。 图5.3.7汽车变速器噪声和振动的功率谱图 利用功率谱估计还可以检查桥桩是否有断裂。在施工中,由于工程质量或地质等原因,桥桩可能断裂,如何发现或判断断裂位置是工程施工中必须解决的问题。当然不可能将桥桩挖出来检查,若用钻机钻出桥桩的混凝土芯来检查,一是费用大,二是有损桥桩强度,应用频谱分析可以很容易解决这个问题。具体方法是,在桩子上端安装一小段水泥管,内部充满水,如图5.3.8所示,在水中插入两个电极,并用仪器对其进行控制。在仪器中的电容充电后,接通线路开关,电极之间产生一个电脉冲。电脉冲通过水,使桩子受到一个垂直于地面的脉冲。同时,在桩子上端装有加速度传感器,测得的加速度信号和处理后得到的频谱如图5.3.9所示。若桥桩没有断裂,则测得的波形比较整齐,只有一个低频的峰值; 若有断裂,则明显地在高频上出现峰值,根据出现的高频成分的频率,很容易估算出断裂的位置,其误差在1m以内。 图5.3.8桥桩试验示意图 图5.3.9桥桩断裂的故障判断 4. 利用实测的荷载谱控制振动台来模拟随机环境 利用实测的荷载谱控制振动台来模拟随机环境。这在研究机件的强度、寿命和可靠性等方面的工作中,得到广泛的应用。 如图5.3.10所示,当汽车在行驶时,将测得的振动加速度信号记录在磁带上,由磁带上的信号回放经分析得到标准的功率谱 P(f),经过逆傅里叶变换和D/A转换,再经过功率放大器,使振动台激振起来; 在振动台上也测出加速度,经过傅里叶变换得到 P^(f); 将 P^(f)与 P(f)进行比较,得到ΔP(f),再来修正功率放大器的信号,使振动台上的振动加速度与实测的汽车振动加速度的功率谱相一致,这样便可以在振动台上模拟道路的随机环境和保证试验工况的稳定性,并可以进行强化快速试验,从而缩短试验时间。 图5.3.10随机环境的模拟系统 5.3.6互谱密度及其估计 设有两个联合平稳随机过程 X(t)和 Y(t),若 x(t,ξ)和 y(t,ξ)分别为 X(t)和 Y(t)的某一个样本函数,相应的截断函数是 xT(t,ξ)和 yT(t,ξ),而 xT(t,ξ)和 yT(t,ξ)的傅里叶变换分别是 XT(Ω,ξ)和 YT(Ω,ξ)。按照与前面对功率谱密度相同的分析方法,定义 X(t)和 Y(t)的互谱密度函数为 PXY(Ω)= limT→∞12TE{XT(-Ω,ξ)YT(Ω,ξ)} (5.3.40) PYX(Ω)= limT→∞12TE{YT(-Ω,ξ)XT(Ω,ξ)} (5.3.41) 对于实函数 xT(t,ξ)和 yT(t,ξ),其频谱一般都是复函数,并有 X*T(Ω,ξ)=XT(-Ω,ξ),Y*T(Ω,ξ)=YT(-Ω,ξ) 所以,上面的两个式子可以写成 PXY(Ω)= limT→∞12TE{X*T(Ω,ξ)YT(Ω,ξ)} (5.3.42) PYX(Ω)= limT→∞12TE{Y*T(Ω,ξ)XT(Ω,ξ)} (5.3.43) 同前,式中 ξ 表示某次试验结果,书写时常略去。由式(5.3.42)和式(5.3.43)可见,互谱密度与实平稳过程的自功率谱密度不同,它不再是Ω的实的、非负的偶函数,而是具有下面所述的性质。 性质1PXY(Ω)=PYX(-Ω)=P*YX(Ω)(5.3.44) 由定义式(5.3.42)和式(5.3.43)即可得证。 性质2Re[PXY(Ω)]和Re[PYX(Ω)]是Ω的偶函数; Im[PXY(Ω)]和Im[PYX(Ω)]是Ω的奇函数。 这可以利用性质1来证明。 性质3若随机过程 X(t)和 Y(t)联合平稳,RXY(τ)绝对可积,则互谱密度和互相关函数构成傅里叶变换对,即 PXY(Ω)=∫∞-∞RXY(τ)e-jΩτdτ(5.3.45) PYX(Ω)=∫∞-∞RYX(τ)e-jΩτdτ(5.3.46) RXY(τ)=12π∫∞-∞PXY(Ω)ejΩτdΩ(5.3.47) RYX(τ)=12π∫∞-∞PYX(Ω)ejΩτdΩ(5.3.48) 以上关系式可以用证明维纳辛钦公式的同样方法证明。 例5.3.4已知平稳过程 X(t)和 Y(t)的互谱密度为 PXY(Ω)=a+jbΩΩ0,-Ω0<Ω<Ω0 0,其他 式中,Ω>0,a,b 为实常数,求互相关函数RXY(τ)。 解 RXY(τ)=12π∫Ω0 -Ω0a+jbΩΩ0ejΩτdΩ =a2π∫Ω0-Ω0ejΩτdΩ+ jb2πΩ0∫Ω0-Ω0ΩejΩτdΩ =1πΩ0τ2{(aΩ0τ-b)sin(Ω0τ)+bΩ0τcos(Ω0τ)} 对有限长离散序列x(n)和 y(m),它们的互谱密度可以写成 Pxy(k)=∑N-1m=0Rxy(m)e-j2πNmk=∑N-1m=01N∑N-1n=0x(n)y(n+m)e-j2πNmk =1N∑N-1m=0x(n)ej2πNnk∑N-1m=0y(m)e-j2πNmk=1NX*(k)Y(k)(5.3.49) 式中,X*(k)是信号x(n)的DFT的共轭; Y(k)是信号 y(m)的DFT; N是有限长离散序列的长度。 互谱密度在信号频谱分析中有很多重要的应用: (1) 利用互谱密度可以得到系统的频率响应函数; (2) 互谱密度能够识别动力学系统的特性; (3) 互谱密度能够确定响应对激励的滞后时间。 5.4谱估计中的几个问题 在对信号进行谱分析过程中,有几个实际的问题需要考虑,下面分别进行介绍。 5.4.1数据预处理 任何信号处理系统总是要对干扰和信号同时进行一番处理,这里所谓信号就是观测者所感兴趣和需要的内容,除此之外都看作干扰或噪声。因此,为了便于实现后续的处理,对采集的数据进行预处理,从而把信号从噪声或干扰中提取出来是非常必要的。 预处理一般还包括对采集数据可疑值的排除,区分周期性函数与随机性数据,确定是否为平稳随机数据等。对于随机数据,如果产生数据的主要物理条件不随时间变化,则可看做是平稳随机信号; 如果信号的功率谱明显地出现尖峰,则认为是周期信号。应根据已知的概率分布规律,对所测得的数据进行检查,从而判断出哪些随机数据是可疑值,然后加以剔除。例如,检验平稳随机过程是否为正态分布,可以将测量数据的概率密度函数与理论正态分布作比较。 预处理的基本方法有预滤波、零均值变换、移动趋势项和信号求平均值等,在实际中,应根据数据的性质和谱估计,选用相应的方法[8]。 1. 预滤波 当信号需要平滑或抑制不需要的频率分量时,可采用滤波的方法。例如,生物信号通过高增益的放大器以后,由于饱和出现低频失真,或数字硬件出现溢出,则可利用高通滤波器进行处理。为了避免因不满足采样定理而出现的频率混叠,可利用低通滤波器来限制原信号的频带宽度,同时还可以减少高频噪声,特别是对于信号求平均的情况。作用于采样信号的数据滤波器,除了可以减少噪声以外,还能抵消漂移和避免当数据通过有限长时窗时所引起的功率泄漏。 2. 零均值变换 为了简化数字信号处理过程的计算,希望通过变换使数据的平均值为0。设原数据为序列x(n),变换后的零均值数据为序列x^(n),则有 x^(n)=x(n)- x-(n)(5.4.1) 式中,x-(n)是原序列x(n)的平均值,即 x-(n)=1N∑Nn=1x(n)(5.4.2) 平均值可以对第1个到第N个的所有采样值进行平均。如果数据是在不同时间段采样得来的,则求平均及零均值变换可在每一段时间内单独进行。 3. 趋势项的移动 被分析的数据中往往包含一个基本趋势,在某些情况下,这个趋势是所需要的,因此,把它分离出来有重要意义。但是,在某些情况下对信号特征的计算会引起误差,因此,在进行预处理时应把它们去掉,特别在进行谱分析的情况下,这样可以避免对峰值的错判及防止对功率谱计算不准确。 趋势可以有几种形式,若趋势由慢变化的信号组成则可以通过数字高通滤波器把它去掉; 如果趋势的周期比记录时间大,则不能采用滤波的方法,在这种情况下,可把它近似看作在整个记录时间内的一个线性变化,采用均方误差最小的方法去掉趋势项。图5.4.1表示去掉趋势项的过程。 图5.4.1去掉趋势项的过程 5.4.2频谱泄漏与窗函数 在进行谱分析时,每次采样所得的时域信号的两端都被截断,如果该信号是确定频率的正弦信号,按理论分析其功率谱,只在确定频率处有一根谱线。但是,实际截取时域信号的是矩形窗函数,结果得到的频谱图有多个频率成分,除了应有的主瓣外,又出现了许多旁瓣。这种谱分析中的失真现象称为泄漏。信号的能量原本应该集中于确定频率处,而现在部分能量泄漏到其他频率上去了。下面对这种现象进行具体分析。 设被处理的信号为x(t)=acos(ω0t),其真实频谱为ω=ω0的一条谱线。但是,对其进行有限时间长度的傅里叶变换时,取 [-T,T]区间代替(-∞,+∞)区间,则有 X(ω)=12π∫T-Tacos(ω0t)e-jωtdt =a2πsin[(ω-ω0)T]ω-ω0+[sin(ω+ω0)T]ω+ω0(5.4.3) 当 T→∞ 时,X(ω)有在 ±ω0处的离散谱线; 当 T 为有限长时,X(ω)呈现出 sinxx 型函数,也就是以 ±ω0为中心,谱线宽度为Δω=2π/T 的主瓣组成的连续谱,原来集中于ω0的功率变成分散在一个与截断长度 T 有关的较宽的频率范围内,但是,总功率不变。 泄漏是由于截取有限时间记录引起的,这使得实际输入信号波形成为非周期性的(如图5.4.2所示),就会产生泄漏[9]。 图5.4.2时间记录中非周期性的输入信号 避免泄漏的最佳方法是,选择合适的窗函数,使信号截断的锐角钝化,使频谱的扩散减到最小。这相当于将时间记录乘上一个函数,这个函数在时间记录的中间最大,在时间记录的两端则为0,从而能将FFT能量集中在时间记录的中间。这类函数称为窗函数。它使我们可以通过一个狭窄的窗口来观察数据。图5.4.3表示窗函数的作用。 时间记录中的非周期性的数据,经过窗口化后,频谱将会明显改善,如图5.4.4所示,其中,图5.4.4(c)采用窗函数的FFT结果显然更接近于正弦波的单条谱线。 为了便于选择,这里给出评价窗函数的4个指标,如图5.4.5所示[4]。 1) 最大旁瓣值比 最大旁瓣值比是最大旁瓣值A旁max与主瓣峰值A峰之比,用对数形式表示。因为最大旁瓣常常使人们分不清它是由于泄漏造成的,还是系统固有的峰值,因此,最大旁瓣比值是评价窗函数的重要指标之一。矩形窗的最大旁瓣比值为-13dB,即 20lg(A旁max/A峰)=-13,A旁max/A峰=0.22 也就是说,最大旁瓣值约为主瓣峰值的 1/5。如果采用后面所示的P310窗,最大旁瓣比值可降低至-93.3dB,即旁瓣值约为主瓣峰值的 2.16×10-5,这时旁瓣值就可忽略不计了。 图5.4.3窗函数及其作用 图5.4.4用窗函数减少泄漏现象 图5.4.5评价窗函数指标的示意图 2) 旁瓣衰减率 旁瓣衰减率以10个相邻旁瓣峰值的衰减比的对数来表示,即为20lg(A旁10/A旁1)。例如,矩形窗或P100哈明窗,其旁瓣衰减率为-20dB,也就是第10个旁瓣峰值为第1个旁瓣峰值的1/10; 而P110汉宁窗的旁瓣峰值为-60dB,即第10个旁瓣峰值为第1个旁瓣峰值的1/1000; P220汉宁窗为-100dB; P330汉宁窗为-140dB。显然,这些窗函数使旁瓣消失很快。 3) 主瓣峰值可能最大误差 ε 如图5.4.6所示,P峰 为理论上计算的连续频谱图上峰值顶点读数,P读为离散化处理时得到的最高峰值谱线处的读数,一般都比 P峰小,则有 ε=P读P峰-1×100%(5.4.4) 图5.4.6主瓣峰值读数误差示意图 由于信号频率是变化的,频率分辨率 Δf 是根据采样频率确定的,而被处理的信号频率很难正好是 Δf 的整数倍,因此,离散化的谱线很难正好在主瓣顶峰处,所以,就产生了误差。例如矩形窗,当其主瓣半个带宽为1Hz时,又设定实际谱线距离顶峰在0.5Hz以内,这时可能的最大误差为-36.34%。 4) 主瓣宽 有时在进行频谱分析时,仅需精确地读出主瓣的频率,这时主瓣窄一些就好分辨,而矩形窗的主瓣最窄,所以,在一些分析中,当对幅值精度要求不高时,可采用矩形窗。 在对实际信号进行分析和处理以及设计数字滤波器时,都会用到不同形式的窗函数。这些窗函数的提出和改进过程,是人们大胆尝试、不断改进的过程[10]。 由于计算时间和计算机资源的限制,我们只能取有限长度的离散信号进行处理。这是一个对信号进行截取的过程,相当于对信号加了个窗函数。若信号原本是无限长的,则只取其中一段; 若原本是周期信号,则不截取整周期就会产生频谱泄漏,造成信号的失真。其原因是: 用窗去截取信号,在时域就是用窗函数与信号相乘,对应的频域就是将窗函数的频谱与原信号的频谱进行卷积。非无限长窗函数的频谱不是脉冲函数,有主瓣和旁瓣,与原信号的频谱卷积后,就会产生频谱泄漏、幅值波动和周期延拓等负面影响。为此,要选取合适的窗函数,以便尽量减小负面影响。 迄今为止,人们已经提出了多种窗函数。最简单就是用矩形窗去截取信号。从频域分析用矩形窗截断信号造成的问题: 矩形窗的频谱是主瓣最窄、旁瓣幅值很大; 主瓣窄有利于计算信号的频率,但是,不利于计算信号的幅值; 旁瓣幅值大造成计算信号幅值时误差很大。从时域分析用矩形窗截断信号带来的问题: 矩形窗的边缘过于陡峭,截断造成原信号周期延拓后,边缘部分的影响所占比重较大。 为此,提出采用三角形窗。与矩形窗相比,这种窗在时域上边缘部分的幅值很小; 其频谱是主瓣宽度变宽,旁瓣幅值较小。 在三角窗的基础上,考虑到傅里叶变换的调制特性,再利用三角函数与指数函数之间的关系,提出了汉宁窗。这种窗函数的频谱可以看成是3个矩形窗频谱的叠加,使得旁瓣互相抵消,能量更加集中于主瓣。所以,与三角形窗的频谱相比,汉宁窗的主瓣宽度相同; 但是,旁瓣小。 再对汉宁窗进行改进,调整了各项之间的比例,提出了哈明窗,其旁瓣幅值更小。 在此基础上,提出布莱克曼窗,增加余弦的二次谐波分量,使得主瓣变宽,旁瓣幅值最小。 这几种窗函数的横向比较和简要分析如表5.4.1所示。分析窗函数的发展和改进过程,我们不难发现,这是人们为了追求窗函数频谱旁瓣更小这个目标,不断进行尝试,探索出来的。 表5.4.1各种窗函数的比较和分析 窗函数时 域 特 点频 谱 特 点原因 矩形窗边缘处幅值很大主瓣窄,旁瓣很大边缘影响大,旁瓣幅值大 三角窗边缘处幅值小瓣宽度变宽,旁瓣幅度较小边缘影响变小,使旁瓣幅值变小 汉宁窗边缘处幅值更小主瓣宽度与上面相同,旁瓣幅度小利用调制特性,使旁瓣互相抵消 哈明窗边缘处幅值较小主瓣宽度与上面相同,旁瓣幅度更小调整各项系数,使旁瓣更小 布莱克曼窗边缘处幅值最小主瓣变宽,旁瓣幅度最小增加余弦二次谐波,使旁瓣最小 5.4.3谱估计的基本步骤 为了保证信号处理的精度和可靠性,对谱估计的实际分析中采用的基本步骤,作如下考虑。 (1) 估计待分析信号的频率范围和频率上限 fm。若需要,可以先对信号进行滤波,去掉过高的频率分量。 (2) 根据分析精度的要求,设定谱分析中的频率分辨率Δf,则有 Δf=1T=1NΔt=fsN(5.4.5) 式中,T 为总采样时间长度; N为采样点数; Δt 为采样时间间隔; fs 为采样频率。 (3) 选定采样间隔 Δt,使采样频率 fs≥2fm。 (4) 由于频域中每个点要由与频率有关的两个值(幅值与相位或实部与虚部)来决定,因此,时域中的N个采样点变换到频域时,只能决定 N/2 个复数点。准确地说,因直流和最高频率只有实数值,故有(N/2)+1个实数点和(N/2)-1个虚数点。谱分析的频带宽Fmax为 Fmax=ΔfN/2(5.4.6) (5) 确定点数 N。当Fmax=fm,则 N=2fm/Δf(5.4.7) 我们希望 Δf 越小越好。由式(5.4.5)可见,当 fs 一定时,N应该大,但是,这会使采样时间长度、存储量和计算量增大。此外,我们常常希望N取2的整次幂。若N已给定,即N不能增加,可用补零的方法使N为2的整次幂。注意,补零不能提高频率分辨率。在进行谱分析时,人们常在有效数据后面补一些零,以达到对频谱做某些改善的目的,这往往会引起一些误解,认为补零会提高分辨率。其理由是,原数据长度为 N1,现数据长度为 N2,由于 Δf1=fs/N1,Δf2=fs/N2, N2>N1 因此, Δf2<Δf1。出现这种错误的原因是把补零后数据的有效长度错认为是 N2,它实际上仍是 N1,即补零不能增加数据的有效长度。虽然理论上补零有一定的优点: ① 可使数据N为2的整次幂,以便于使用FFT算法; ② 补零起到对原DFT的 X(k)作插值的作用,一方面,克服“栅栏”效应,使谱的外观得到平滑; 另一方面,由于对数据截断时所引起的频谱泄漏,有可能在频谱中出现一些难以确认的峰值,补零后有可能消除这种现象。 但是,在实际中应该尽量多采集一些数据,避免补零。 5.4.4频谱校正方法 采用周期图谱分析方法进行信号处理,在科研、工程和生产中有十分广泛的应用。但是,这种方法存在以下局限性。 (1) 计算机只能对有限样本进行处理,FFT谱分析也只能在有限区间内进行。由于时域截断产生的能量泄漏,造成谱峰值变小,精度降低。 (2) 采样频率不可能是信号频率的整数倍,而FFT的频谱是离散的,若信号频率在两条谱线之间,这时由峰值谱线反映的频率、幅值和相位就存在较大误差。例如,当加矩形窗且进行非整周期采样时,频谱分析的幅值误差最大达 36.4%,频率误差最大为 0.5Δf,Δf 为频率分辨率。 为此,人们提出了一些校正频谱分析中误差的方法,以满足实际应用对频谱分析精度的要求。下面仅介绍频谱重心校正法,并给出其处理步骤、具体算法和特点[7,11]。 1) 基本原理 频谱重心校正法的基本原理是,利用窗函数主瓣内的谱线求主瓣中心的坐标,得到准确的频率、幅值和相位,根据主瓣函数的特点用重心法规求中心坐标。不同的窗函数有不同的主瓣形状,下面仅以矩形窗为例。 2) 矩形窗的修正公式 矩形窗为 w(k)=1,k=0,1,…,N-1(5.4.8) 其离散傅里叶变换为 W(ω)=∑N-1k=0e-j2πkn/N=1-cos(Nω)+jsin(Nω)1-cosω+jsinω =sin(Nω/2)sin(ω/2)e-jω(N-1)/2(5.4.9) 式中,ω=2πk/N。 矩形谱的模函数为 W0(n)=sin(πn)sin(πn/N)(5.4.10) 当N1时,1/N→0,sin(πn/N)≈nπ/N,故在主瓣区间有 W0(n)=Nsin(πn)πn(5.4.11) 当ω=2πn/N(n=±1,±2,…)时,W=0,所以,主瓣宽度为 4π/N。又因为谱线间隔为 2π/N,所以,主瓣内有两条谱线,分别为第 n 条和第 n+1 条。利用式(5.4.11),并用 Y 代替 W0,则有 Y(n)=sin(πn)πn(5.4.12) 因为 nY(n)+(n+1)Y(n+1)=nNsin(πn)/(nπ)+(n+1)Nsin(nπ+π)/[(n+1)π] =Nsin(πn)-Nsin(πn)=0(5.4.13) 所以,式(5.4.13)说明两条谱线的重心为主瓣中心,按重心法规求中心坐标 x0,则有 x0=nY(n)+(n+1)Y(n+1)Y(n)+Y(n+1)=n+Y(n+1)Y(n)+Y(n+1)(5.4.14) 令 Δn=Y(n+1)Y(n)+Y(n+1)(5.4.15) 由频率的一般形式 f=nfs/N,得到修正公式为 f=(n+Δn)fs/N(5.4.16) 同理,可以推出用主瓣内相邻谱峰最高的两条谱线(第n-1条和第n条)校正频率的公式。 设主瓣峰值为A,则 Y(n)=Asin[π(n-x0)]π(n-x0)(5.4.17) 式(5.4.17)相当于式(5.4.12)乘以系数 A,并平移到 x0处。x0和 A 分别是分析信号的频率和幅值,Y(n)和 Y(n+1)是主瓣内的两条谱线。 将 x0,Y(n)代入式(5.4.17),得到幅值校正公式为 A=Y(n)Δnπsin(Δnπ)(5.4.18) 由式(5.4.9)可知,矩形窗频谱函数的相位角为 φ=-N-12ω(5.4.19) 将ω=2πn/N代入式(5.4.19),在主瓣内近似认为(N-1)/N≈1,则 φ= -nπ,校正量 Δφ=-Δnπ,复频信号与窗频谱函数进行卷积运算时,相位角相加,得到相位校正公式为 θ=arctan(In/Rn)-Δnπ(5.4.20) 式中,Rn 和 In 分别为FFT的实部和虚部。 3) 几点讨论 ① 该校正方法计算速度快,精确度较高; ② 对于复杂的窗函数校正公式的推导困难; ③ 当频率间隔过小、主瓣重叠时,该方法不适用。 5.5平稳随机信号通过线性系统 设 X(n)为一平稳随机信号,当它通过一个线性移不变系统 H(z)后,输出为Y(n),即 Y(n)=X(n)*h(n)=∑∞k=-∞ X(k)h(n-k) 所以,Y(n)也是随机的,且也是平稳的。若 X(n)是确定性信号,则 Y(ejω)=X(ejω)H(ejω) 由于随机信号不存在傅里叶变换,因此,需要从相关函数和功率谱的角度来研究随机信号通过线性系统的响应。为了讨论问题方便起见,假设 X(n)是实信号,这样,Y(n)也是实的,输入 X(n)和输出 Y(n)有以下关系。 (1) RY(m)=RX(m)*h(m)*h(-m)(5.5.1) 因为 RY(m)=E{Y(n)Y(n+m)} =E∑∞k=-∞h(k)X(n-k)∑∞i=-∞h(i)X(n+m-i) =∑kh(k)∑ih(i)E{X(n-k)X(n+m-i)} =∑kh(k)∑ih(i)RX(m+k-i)=∑kh(k)R(m+k) 式中, R(m+k)=h(m+k)*RX(m+k) 令 m+k=l,则 RY(m)=∑∞k=-∞h(k)R(m+k)=∑∞l=-∞R(l)h(l-m)=R(m)*h(-m) 所以, RY(m)=RX(m)*h(m)*h(-m) (2) PY(ejω)=PX(ejω)H(ejω)2(5.5.2) 对式(5.5.1)两边作傅里叶变换,就可得到式(5.5.2)。 (3)RXY(m)=RX(m)*h(m)(5.5.3) RXY(m)=E{X(n)Y(n+m)}=EX(n)∑∞k=-∞h(k)X(n+m-k) =∑kh(k)E{X(n)X(n+m-k)} =∑kh(k)RX(m-k)=RX(m)*h(m) (4)PXY(ejω)=PX(ejω)H(ejω)(5.5.4) 对式(5.5.3)两边取傅里叶变换,就得到式(5.5.4)。 例5.5.1一个简单的两点差分器可以用下式来描述: Y(n)=12[X(n)-X(n-2)] 它可以用来近似计算信号的斜率。设 X(n)为一均值为0、方差为σ2 的白噪声,试求输出 Y(n)的自相关函数和功率谱。 解因为 X(n)为一白噪声序列,所以 RX(m)=σ2Xδ(m),PX(ejω)=σ2X 由所给系统得 H(ejω)=je-jωsinω=e-j(ω-π/2)sinω h(0)=1/2,h(2)=-1/2,h(n)=0,n为除0,2以外的其他值。 由式(5.5.1)可得 RY(m)=σ2X/2,m=0 -σ2X/4,m=±2 0,其他 以上结果也可以直接由RY(m)的定义求出。对RY(m)求傅里叶变换,有 PY(ejω)=σ2Xsin2ω 5.6与本章内容有关的MATLAB函数 与本章内容有关的MATLAB函数主要用于功率谱估计、频率响应分析和相关函数计算。 1. periodogram() 该函数采用周期图估计一个信号的功率谱密度或功率谱,直接将信号的采样序列进行傅里叶变换求取功率谱估计,其调用格式为 [Pxx,f]=periodogram(xn,nfft,fs,window) [Pxx,f]=periodogram(xn,nfft,fs,window,'range') 式中,xn是信号序列; nfft是fft的长度; fs是采样频率; window是选用的窗函数,窗长要小于或等于nfft; 输出的Pxx是估计出的功率谱,f是频率坐标。默认值nfft=256,fs=1Hz,采用矩形窗。range是频率范围选项: twoside是计算双边功率谱,oneside是计算单边功率谱。periodogram()函数允许在参数中指定一个空矩阵,即[ ],表示使用函数的默认参数设置。 2. pwelch() 该函数是用Welch法估计一个信号的功率谱,Welch法就是在重叠的加窗分段信号上作平均周期图,其调用格式为 [Pxx,f]=pwelch(xn,nfft,fs,window,noverlap) [Pxx,f]=pwelch(xn,nfft,fs,window,noverlap,'range') Pxx=pwelch(xn,window,noverlap) 式中,xn是信号序列; nfft是fft的长度; fs是采样频率; window是选用的窗函数,窗长要小于或等于nfft; noverlap是估计功率谱时每一段重叠的长度,要小于nfft; 输出的Pxx是估计出的功率谱; f是频率坐标。默认值nfft=256,fs=1Hz,50%重叠(即前一段信号和后一段信号有一半是重叠的),采用汉宁窗。range的意义同periodogram()函数。pwelch()函数也允许在参数中指定一个空矩阵,即[ ],表示使用函数的默认参数设置,例如, [Pxx,f]=pwelch(xn,[ ],fs,[ ],noverlap) 表示nfft=256,采用汉宁窗。 3. psd() 该函数是用Welch法估计一个信号的自功率谱,其调用格式为 Pxx=psd(xn) Pxx=psd(xn,nfft) [Pxx,f]=psd(xn,nfft,fs) Pxx=psd(xn,nfft,fs,window) Pxx=psd(xn,nfft,fs,window,noverlap) Pxx=psd(xn,nfft,fs,window,noverlap,'dflag') 式中,Pxx为信号的自功率谱密度估计; dflag 为去除信号趋势量的选项: line表示去除直线趋势量,mean表示去除均值量,none表示不作去除信号趋势处理; xn,nfft,fs,window,noverlap的含义同pwelch()函数。默认值nfft=256,fs=2,noverlap=0,采用汉宁窗。若xn为实序列,则Pxx只计算频率为正的功率谱估计; 若xn为复序列,则Pxx计算正负频率的功率谱估计。psd()函数允许在参数中指定一个空矩阵[ ],表示使用函数的默认参数设置。 例5.6.1计算叠加白噪声序列的正弦信号 2sin(0.12πn)+sin(0.28πn)的功率谱。 解MATLAB程序如下: %examp5_1.m n=0:1000; xn=2*sin(0.12*pi*n)+sin(0.28*pi*n)+randn(size(n)); %产生信号 nfft=input('the size of fft='); %输入fft运算长度 window=hamming(256); noverlap=input('the amount of overlap=');%输入覆盖点数 [Pxx,f]=psd(xn,nfft,2,window,noverlap);%计算功率谱密度 plot(f/2,10*log10(Pxx));grid%绘图 xlabel('\omega/\pi'); %加标注 ylabel('Power Spectrum,dB'); %nfft=1000,noverlap=128 结果如图5.6.1所示。 图5.6.1功率谱估计 ■ 4. spectrum() 该函数也是用Welch法估计一个信号的功率谱,调用格式为 [Pxx,f]=spectrum(xn,nfft,fs,window,noverlap,'dflag') 式中,参数含义同psd()函数,使用方法也同psd()函数。 5. csd() 该函数是用Welch法估计两个等长信号x和y的互功率谱,csd()函数用信号x和y的FFT乘积来形成周期图,所以其结果是复数。其调用格式为 Pxy=csd(x,y) Pxy=csd(x,y,nfft) [Pxy,f]=psd(x,y,nfft,fs) Pxy=csd(x,y,nfft,fs,window) Pxy=csd(x,y,nfft,fs,window,noverlap) Pxy=csd(x,y,nfft,fs,window,noverlap,'dflag') 式中,x,y 为两个信号序列; Pxy为x,y的互功率谱估计; 其他参数的意义同psd函数。 小结 本章介绍了随机信号分析与处理的有关内容,着重讲述了随机信号的相关函数分析及应用,重点讨论了功率谱估计,其中包括功率谱密度的定义、性质、与自相关函数的关系、估计方法和应用,简要介绍了互谱分析。本章还阐述了谱估计中的几个重要问题,包括预处理、频谱泄漏与窗函数之间的关系、谱估计的基本步骤及频谱校正方法。最后,本章介绍了随机信号通过线性移不变系统的响应。 习题和上机练习 5.1功率谱在实际中能解决什么问题? 5.2若随机过程 X(t)的自相关函数为 RX(τ)=1-τ,τ≤1 0,τ>1 求功率谱密度。 5.3求 Y(t)=X(t)cos(Ω0t+Φ)的功率谱密度。式中,X(t)为平稳随机过程,Φ为(0,2π)上均匀分布的随机变量,Ω0为常数,X(t)与Φ互相独立。 5.4已知平稳过程 X(t)的自相关函数为 RX(τ)=1-τT,-T≤τ<T 0,其他 求 X(t)的功率谱密度 PX(Ω),并画出其图形。 5.5推导用周期图法求互功率谱密度估值的过程,画出框图。 5.6某平稳随机信号最高频率为50Hz,现以200Hz进行采样,若要求频率分辨率为2Hz,试确定进行谱估计时每段数据的长度。 5.7已知平稳随机信号的一个观测序列x(n)=[1,1,0,0],试用周期图法求出其功率谱的估值。 5.8设信号频率为1kHz,采样频率为10kHz,采样数据个数为1024,试求谱估计的频率分辨率和可分析的频率范围; 若要提高频率分辨率,应采取什么措施? 5.9在经典谱估计中利用直接法进行功率谱估计时,补零起什么作用?为什么说补零不能提高频率分辨率? 5.10试对频率混叠和频谱泄漏进行讨论,指出它们产生的原因,并给出解决方法。 5.11一段记录包含N点采样,其采样频率 fs=1000Hz。用平均法改进周期图估计时将数据分成互不交叠的K段的方式,其中每段数据长度M=N/K。假定在频谱中有两个相距为0.04π(rad)的谱峰,要分辨它们,M应取多大? 5.12在用某台FFT仪做谱分析时,选用抽样点数N必须是2的整数次幂。已知待分析的信号中,上限频率≤1.25kHz。要求谱分辨率≤5Hz。试确定以下参数: (1) 一个记录中的最少抽样点数; (2) 相邻样点间的最大时间间隔; (3) 信号的最小记录时间。 图题5.13 5.13若系统的输入过程 X(t)为平稳过程,系统的输出(见图题5.13)为 Y(t)=X(t)+X(t-τ) 求证 Y(t)的功率谱密度为 PY(Ω)=2PX(Ω)[1+cos(Ωτ)] 5.14假设有一个二阶AR系统的系统函数是 H(z)=11-2ρz-1cosθ+ρ2z-2,ρ<1 当一个方差为σ2w 的白噪声序列 w(n)输入该系统中时,求输出序列 y(n)的频谱,并分析 ρ 和 θ 的变化对输出频谱的影响。 *5.15试求下列差分方程所描述的输出序列x(n)的功率谱并作图: (1) x(n)=-0.81x(n-2)+w(n)-w(n-1) (2) x(n)=w(n)-w(n-2) (3) x(n)=-0.81x(n-2)+w(n) 式中,w(n)是方差为σ2w(如σ2w=1/12)的白噪声。 *5.16一序列x(n)是由两个频率相距为 Δf 的模拟信号采样得来的,即 x(n)=sin(2π×0.135×n)+cos[2π(0.135+Δf)n]n=0,1,…,15 已知序列长度N=16,试采用周期图法,应用DFT分别计算当 Δf=0.06及 Δf=0.01时的功率谱估计,并通过作图说明从功率谱估计的分布是否能分辨出这两个正弦信号的真实频谱。若N=64又有什么变化? *5.17用MATLAB产生256点的白噪声序列,应用Welch法估计其功率谱,每段长64点,重叠32点,作出输出平均后的功率谱图以及对256点一次求周期图的功率谱图。 *5.18离散信号序列x(n)=2sin(2πfn/fs),长度N=512,fs=32Hz,令 f 取值分别为16Hz和15Hz,计算序列的功率谱,比较两个功率谱图的差别。试采用不同的窗函数,再比较功率谱图的变化。 *5.19已知一个被白噪声 r(t)污染的信号x(t)为 x(t)=2sin(2πf1t)+0.5sin(2πf2t)+0.5sin(2πf3t)+r(t) 其中,f1=25Hz,f2=75Hz,f3=150Hz。应用Welch法进行功率谱估计并绘制功率谱图。 参考文献 [1]应怀樵.波形和频谱分析与随机数据处理[M].北京: 中国铁道出版社,1983. [2]王永德,王军.随机信号分析基础[M].2版.北京: 电子工业出版社,2003. [3]胡广书.数字信号处理理论、算法与实现[M].2版.北京: 清华大学出版社,2003. [4]黄世霖.工程信号处理[M].北京: 人民交通出版社,1986. [5]舒张平,徐科军.互相关分析中有偏估计与无偏估计的比较[J].电气电子教学学报,2017,39(3): 3234. [6]宗孔德,胡广书.数字信号处理[M].北京: 清华大学出版社,1988. [7]徐科军,全书海,王建华.信号处理技术[M].武汉: 武汉理工大学出版社,2001. [8]吴湘淇,聂涛.数字信号处理技术及应用[M].北京: 中国铁道出版社,1986. [9]沈兰荪.高速数据采集系统的原理与应用[M].北京: 人民邮电出版社,1995. [10]徐科军.以能力为导向讲授研究生信号处理课[J].电气电子教学学报,2019,41(2): 5659,76. [11]谢明,丁康. 频谱分析的校正方法[J]. 振动工程学报,1994,7(2): 172178.