1、2.2.7.3 相关分析一、相关 二、互相关函数与自相关函数 三、相关函数的工程意义及应用 一、相关(correlation)相关：用来描述一个随机过程自身在不同时刻的状态间，或者两个随机过程在某个时刻状态间线性依从关系的数字特征。图2.52 变量x和y的相关性 （a）精确相关 （b）中等程度相关 （c）不相关 评价变量x和y间线性相关程度的经典方法：n协方差xy：式中，E表示数学期望值；x=Ex为随机变量x的均值；y=Ey为随机变量y的均值；n相关函数xy：式中x、y分别为x、y的标准偏差，而x和y的方差x2和y2则分别为 NiyixiNyxxyyxNyxE11lim（2.142）11xyy
2、xxyxy（2.143）22xxxE22yyxE（2.144）（2.145）利用柯西许瓦兹不等式(Cauchy-Schwarz inequality)可知|xy|1。当xy=1时，所有数据点均落在y-y=m(x-x)的直线上，因此x,y两变量是理想的线性相关。当xy=0时，(xi-x)与(yi-y)的正积之和等于其负积之和，因而其平均积xy为0，表示x,y之间完全不相关。222yxyxyExEyxE（2.146）二、互相关函数与自相关函数 对于各态历经过程，可定义时间变量x(t)和y(t)的互协方差(cross-covariance)函数为 式中称x(t)与y(t)的互相关(cross-cor
3、relation)函数，自变量称为时移。yxxyTxxTyxxyRdttytxTtytxEC01lim（2.147）TTxydttytxTR01lim（2.148）当y(t)x(t)时，得自协方差(auto-covariance)函数 其中 称为x(t)的自相关(auto-correlation)函数。周期函数的自相关函数仍为周期函数，且两者的频率相同，但丢掉了相角信息。同频相关，不同频不相关。201limxxTxxTxRdttxtxTC（2.149）TTxdttxtxTR01lim（2.150）图2.53 典型的自相关函数和互相关函数曲线(a)自相关函数 （b）互相关函数 例1 求正弦函数x
4、(t)=Asin(t+)的自相关函数。解：正弦函数x(t)是一个均值为零的各态历经随机过程，其各种平均值可用一个周期内的平均值来表示。令t+=,则dt=d/,由此得正弦函数的自相关函数是一个与原函数具有相同频率的余弦函数，它保留了原信号的幅值和频率信息，但失去了原信号的相位信息。自相关函数可用来检测淹没在随机信号中的周期分量。0000)(sin)sin(1)()(1lim)(TTTxdtttATdttxtxTRcos2)sin(sin2)(2202AdARx自相关和互相关函数的估计 和 具有限个数据点N的相关函数估计的数字处理表达式则为：)(xR)(xyR TxdttxtxTR01 Txydt
5、tytxTR01（2.160）（2.161）11NonxrnxnxNrR 11NonxyrnynxNrRNrr,2,1,0（2.162）（2.163）三、相关函数的工程意义及应用 不同类别信号的辨识 图2.55 典型信号的自相关函数 相关滤波(filtering by correlation)图4.79 相关滤波频谱分析仪原理框图 相关测速和测距 图2.56 相关法测量声传播距离 图2.57 带钢测速系统 测量流速和流量 图2.58 相在法测定流量 2.2.7.4 功率谱分析 2.2.7.4 功率谱(power spectrum)分析 一、自功率谱密度函数 二、巴塞伐尔（Parseval）定理
6、 三、互功率谱密度函数 四、自谱和互谱的估计五、工程应用 一、自功率谱密度函数 设x(t)为一零均值的随机过程，且x(t)中无周期性分量，则其自相关函数Rx()在当时有 该自相关函数Rx()满足傅里叶变换的条件 。对它作傅里叶变换可得 其逆变换为 0)(xRdRx)(deRfSfjxx2（2.167）dfefSRfjxx2（2.168）Sx(f)为x(t)的自功率谱密度函数(auto power density spectrum)，简称自谱或功率谱。功率谱Sx(f)与自相关函数Rx()之间是傅里叶变换对的关系，亦即式（2.167）和（2.168）称为维纳辛钦（Wiener-Khintchine
7、）公式。由于Rx()为实偶函数，因此亦为Sx(f)实偶函数。)()(fSRxFTIFTx图2.59 单边功率谱和双边功率谱 当=0时，根据自相关函数Rx()和自功率谱密度函数Sx(f)的定义，可得Sx(f)曲线下面和频率轴所包围的面积即为信号的平均功率；Sx(f)就是信号的功率谱密度沿频率轴的分布，故也称为功率谱。dffSdttxTRxTTTx2221lim0（2.169）二、巴塞伐尔（Parseval）定理设有变换对：按频域卷积定理有 令k=0，有 又令h(t)=x(t)，得)()(fXtx)()(fHth)(*)()()(fHfXthtxdffkHfXdtethtxktj)()()()(2
8、dffHfXdtthtx)()()()(dffXfXdttx)()()(2x(t)为实函数，故X(-f)=X*(f)，于是有v巴塞伐尔定理：信号在时域中计算的总巴塞伐尔定理：信号在时域中计算的总能量等于它在频域中计算的总能量。能量等于它在频域中计算的总能量。式（2.170）称信号能量等式。|X(f)|2称能量谱，是沿频率轴的能量分布密度。在整个时间轴上信号的平均功率可计算为自谱密度函数与幅值谱之间的关系:(rf.(2-169)dffXdffXfXdttx2*2（2.170）dffXTdttxTPTTTT22221lim1lim（2.171）21limfXTfSTx（2.172）对于单边(one
9、-sided)功率谱G(f)也应满足巴塞伐尔定理，故有由此规定 Gx(f)的图形如图2.59中所示。0dffGdffSPxx（2.173）0)(2)(ffSfGxx图2.59 单边功率谱和双边功率谱 根据信号功率（或能量）在频域中的分布情况，将随机过程区分为窄带随机、宽带随机和白噪声等几种类型。窄带过程的功率谱（或能量）集中于某一中心频率附近，宽带过程的能量则分布在较宽的频率上，而白噪声过程的能量在所分析的频域内呈均匀分布状态。三、互功率谱密度函数 若互相关函数Rxy()满足傅里叶变换的条件 ，则定义Rxy()的傅里叶变换 为信号x(t)和y(t)的互功率谱密度函数，简称互谱密度函数(cros
10、s power spectrum)或互谱。根据维纳辛钦关系，互谱与互相关函数也是一个傅里叶变换对，即 因此Sxy(f)的傅里叶逆变换为：dRxy)(dfeRSftjxyxy2（2.175）)()(fSRxyFTIFTxy dfefSRfjxyxy2（2.176）定义信号x(t)和y(t)的互功率为因此互谱和幅值谱的关系为 正如Ryx()Rxy()一样，当x和y的顺序调换时，Syx()Sxy()。但根据Rxy(-)=Ryx()及维纳辛钦关系式，不难证明：其中 dffXfYTdttytxTPTTTT*221lim1lim（2.177）fXfYTfSTxy*1lim（2.178）fSfSfSyxxy
11、xy*（2.179）)()(1lim*fYfXTfStxy Sxy(f)也是含正、负频率的双边互谱，实用中也常取只含非负频率的单边互谱Gxy(f)，由此规定 自谱是f的实函数，而互谱则为f的复函数，实部Cxy(f)称为共谱(cospectrum)，虚部Qxy(f)称为重谱(quad spectrum)，即写为幅频和相频的形式：02ffSfGxyxy（2.180）fjQfCfGxyxyxy（2.181）fCfQarctgffQfCfGefGfGxyxyxyxyxyxyfjxyxyxy22（2.182）四、自谱和互谱的估计定义功率谱亦即自谱的估计值 互谱的估计为 21fXTfSx（2.183）fY
12、fXTfSxy*1（2.184）fXfYTfSyx*1（2.185）五、工程应用 1.求取系统的频响(frequency response)函数 线性系统的传递函数H(s)或频响函数H(j)十分重要，在机器故障诊断等多个领域常要用到它。w例1：机器由于其轴承的缺陷而在机器运行中会造成冲击脉冲信号，此时若用安装在机壳外部的加速度传感器来接收时，必须考虑机壳的传递函数。w例2：当信号经过一个复杂系统被传输时，系统各环节的传递函数便必须要加以考虑。一个线性系统的输出y(t)等于其输入x(t)和系统的脉冲响应h(t)的卷积，即 根据卷积定理，上式在频域中化为式中H(f)即为系统的频响函数。thtxty
13、*（2.192）fXfHfY（2.193）通过自谱和互谱来求取H(f)：对式（2.193）两端乘以各自的复共轭并取期望值有 上式反映出输入与输出的功率谱密度和频响函数间的关系；式中没有频响函数的相位信息，因此不可能得到系统的相频特性。fSfHfSxy2（2.194）如果在式（2.193）两端乘以x(f)的复共轭并取期望值，则有 由于Sx(f)为实偶函数，因此频响函数的相位变化完全取决于互谱密度函数的相位变化。式（2.195）将输入、输出的相位关系完全保留了下来，且在这里输入的形式并不一定限制为确定性信号，也可以是随机信号。fXfXfHfXfY*fSfHfSxxy（2.195）通常一个测试系统往
14、往受到内部和外部噪声的干扰。从而输出也会带入干扰。输入信号与噪声是独立无关的，因此它们的互相关为零。结论：在用结论：在用和和求取系统频响函求取系统频响函数时不会受到系统干扰的影响。数时不会受到系统干扰的影响。(优点优点)2.旋转机械振动特性检测 n旋转机械的转轴部件从起动、升速到额定转速的过程共经历了全部转速的变化，因此在各个转速下的振动状态可用来对机器的临界转速、固有频率和阻尼比等各参数进行辨识。n起动和停车过程则包含了丰富的信息。是常规运行状态下所无法获得的。n“瀑布图(waterfall plot)法”：在机械振动或停车过程中将不同转速下振动的功率谱图迭加而形成的一种图。图2.61 旋转
15、机械的瀑布图 由图可见机器的回转频率n(r/min)及其各次谐波下谱峰高度，由此来得出机器的临界转速、固有频率及阻尼比等数据。从图可见，机器临界转速约为4000r/min，机器振动的高次谐波分量很小，主要是回转频率处的谱峰，因此可判断转子存在有较严重的失衡。此外还可看到图中频率60HZ处有一谱峰值，它不随转速升高而改变，判断为电源的脉动干扰。2.3 数字信号处理数字信号处理(digital signal processing)：利用计算机或专用信号处理设备，以数值计算的方法对信号作采集、变换、综合、估值与识别等处理。一、离散傅里叶变换（DFT）二、离散傅里叶变换的性质 三、采样定理 四、泄漏与
16、加窗处理 五、栅栏效应六、快速傅里叶变换（FFT）一、离散傅里叶变换（DFT）对于一个非周期的连续时间信号x(t)来说，它的傅里叶变换应该是一个连续的频谱X(f)，其运算公式根据第二章的内容有 dtetxfXFTftj2:dfefXtxIFTftj2:（2.199）（2.200）图2.63 傅里叶变换的几种类型 对于无限连续信号的傅里叶变换共有四种情况：n对于非周期连续信号X(t)，频谱X(f)是连续谱；n对于周期连续信号，傅里叶变换转变为傅里叶级数，因而其频谱是离散的；n对于非周期离散信号，其傅里叶变换是一个周期性的连续频谱；n对于周期离散的时间序列，其频谱也是周期离散的。结论：若x(t)是
17、周期的，频域中X(f)必然是离散的，反之亦然。若x(t)是非周期的，则X(f)一定是连续的，反之亦然。第四种亦即时域和频域都是离散的信号，且都是周期的，给我们利用计算机实施频谱分析提供了一种可能性。对这种信号的傅里叶变换，我们只需取其时域上一个周期（N个采样点）和频域一个周期（同样为N个采样点）进行分析，便可了解该信号的全部过程。DFT的定义：对有限长度的离散时域或频域信号序列进行傅里叶变换或逆变换，得到同样为有限长度的离散频域或时域信号序列的方法，便称为离散傅里叶变换（discrete Fourier transform,DFT）或其逆变换（IDFT）。离散傅里叶变换的公式：式中x(n)和X
18、(k)分别为 和 的一个周期，此处将t和f0均归一化为1。nkNNonNonnkNjWnxenxkX112 11211NoKNoKnkNnkNjWkXNekXNnx（2.205）（2.206）NjNeW2)(tnx)(0kfX离散傅里叶变换意义：可以对任意连续的时域信号进行采样和截断并对其作离散傅里叶变换的运算，得到离散的频谱，该频谱的包络即是对原连续信号真正频谱的估计。离散傅里叶变换的过程：n时域采样(sampling in t-domain)；n时域截断(truncation in t-doman)；n频域采样(sampling in f-domain)。图2.64 离散傅里叶变换的图解过
19、程（一）图2.64 离散傅里叶变换的图解过程（二）图2.64 离散傅里叶变换的图解过程（三）1.线性性线性性如果则式中 a 和 b 为常数.2.3.2 离散傅立叶变换的性质离散傅立叶变换的性质)()()()(2211kXnxDFTkXnxDFT kbXkaXnbxnaxDFT2121(2.214)2.序列的移动性序列的移动性 如果时移:频移:2.3.2 离散傅立叶变换的性质离散傅立叶变换的性质)()(kXnxDFT NkmjekXmnxDFT/2(2.215)mkXenxDFTNkmj/2(2.216)3.对称性对称性1)如果 x(n)为复序列,且则2)如果 x(n)为实序列,且则式中 XR(
20、k)，XI(k)分别为X(k)的实部和虚部.2.3.2 离散傅立叶变换的性质离散傅立叶变换的性质)()(kXnxDFT kXnxDFT*(2.217)()(kXnxDFT kXkXkNXkXkNXkXkXkNXkXkXkNXkXkXIIIRRRargarg*(2.218)3)如果x(n)为偶序列,即,x(n)=x(-n),则 X(k)为偶序列.4)如果x(n)为奇序列,即,x(n)=-x(-n),则 X(k)为纯虚序列.2.3.2 离散傅立叶变换的性质离散傅立叶变换的性质4.巴什伐尔定理巴什伐尔定理5.卷积卷积若则2.3.2 离散傅立叶变换的性质离散傅立叶变换的性质 11221NonNokkX
21、Nnx(2.219)()()()(kYnyDFTkXnxDFT kYkXnynx*(2.220)kYkXnynx*(2.221)三、采样定理(sampling theorem)混叠（aliasing）：若采样率过低即采样间隔大，则系列的离散时间序列可能不能真正反映原始信号的波形特征，在频域处理时会出现频率混淆。图2.65 不同采样率对采样信号产生的影响（一）图2.65 不同采样率对采样信号产生的影响（二）采样定理(sampling theorem)：为避免混叠产生，要求的采样频率fs必须高于信号频率成分中最高频率fmax的两倍，即 乃奎斯特（Nyquist）频率：在给定的采样频率fs条件下，信
22、号中能被分辨的最高频率。只有低于乃奎斯特频率的频率成分才能被精确地采样，亦即为避免频率混淆，应使被分析信号的最高频率fmax低于乃奎斯特频率。max2ffs（2.222）2sNyqff（2.223）图2.66 混叠产生的条件 四、泄漏(leakage)与加窗(windowing)图2.67 余弦信号加窗截断造成的泄漏现象 抑制或减小泄漏效应的方法：选择性能更好的特殊窗来替代矩形窗，亦即加窗处理。评价窗函数的性能指标：1.3dB带宽B：它是主瓣归一化的幅值下降至-3dB时的带宽。归一化|W(f)|=20lg|W(f)/W(0)|，带宽B的单位为或f。2.旁瓣幅度A(dB)，表示为最大旁瓣峰值As
23、max与主瓣峰值Am之比，即20lg(Asmax/Am)。3.旁瓣峰值衰减率D（dB/decade），表示为最大旁瓣峰值与相距十倍频处的旁瓣峰值之比，也是以分贝表示。理想的窗函数应具有最小的B和A以及最大的D。（图2.67）图2.69 常用窗函数的时域图像 图2.70 常用窗函数的频谱 五、栅栏效应(picket fence effect)栅栏效应：若信号中某频率成分的频率fi等于k/T，即它与输出的频率采样点相重合，那么该谱线便可被精确地显示出来；反之若fi与频率采样点不重合，便得不到显示，所得的频谱便会产生误差。频率分辨率f：两条谱线间的距离。当被分析的时域信号长度T（即窗宽T=NTs）和
24、采样频率fs被确定之后，则频率分辨f也被确定：当N增大，即频率分辨率f增大后，是否一定得到精确的离散频谱呢？TNffs1（2.231）例：对余弦信号cos2f0t作DFT。图2.72 周期信号作整周期截取的DFT（一）图2.72 周期信号作整周期截取的DFT（二）图2.73 周期函数作非整周期截取的DFT 结论：对周期信号作整周期截取是获取正确频谱的先决条件。六、快速傅里叶变换（FFT）离散傅里叶变换的计算公式为：式中 N个点的X(k)需做N2次复数乘法和N(N-1)次复数加法。而做一次复数乘法需要做四次实数相乘和两次实数相加，做一次复数加法需要做两次实数相加。例：N=1024时，则需要总共1
25、,048,576次复数乘，即4,194,304次实数乘法。112:NonnkNNonnkNjWnxenxkXDFT 11211:NoknkNNoknkNjWkXNekXNnxIDFT1,1,0Nk1,1,0Nn（2.232）（2.233）NjNeW2快速傅里叶变换(FFT，Fast Fourier Transform)算法的本质:充分利用因子WN的周期性和对称性。n对称性：n周期性：因子WN中的N2个元素只有个N个独立值，其中N/2个值与其余N/2个值数值相等，符号相反。FFT算法的基本思想：避免运算中的重复运算，将长序列的DFT分割为短序列的DFT的线性组合，从而达到整体降低运算量的目的。效
26、果：使原来的N点DFT的乘法计算量从N2次降至为N/2log2N次，如N=1024，则计算量现在为5120次，仅为原计算量的4.88%。nkNNnkNWW2nkNnkNNWW（2.234）（2.235）时间抽取(decimation-in-time)基2算法 对式(2.232)，令N=2M,将x(n)序列分割成长度各为N/2的奇序列和偶序列，即令n=2r和n=2r+1,，r=0,1,N/2-1则式（2.232）重写为式中 这是因为 1221221212122122122NorrkNNorkNrkNNorNorrNrkNWrxWWrxkWrxWrxkX（2.236）NjNjNeeW4222222
27、222NNjNjNWeeW令则式（2.236）可改写为而 因此将式（2.239）完整地写成 12212,1,02NorrkNNkWrxkA 12212,1,012NorrkNNkWrxkB（2.237）（2.238）12/,1,0NkkBWkAkXkN（2.239）kBkNBkAkNA22 12,2,1,022NkkBWkANkXkBWkAkXNkNkN（2.240）又因为 ，因此最终可得 总计为N(N+1)/2次复数乘法，工作量减少一半。12NNW 12,2,1,02NkkBWkANkXkBWkAkXkNkN（2.241）图2.74 分割一次后的A(k)、B(k)及X(k)之间的关系（N=8
28、）按照上述思路继续对A(k)和B(k)作奇偶序列分解。令r=2l，r=2l+1，l=0,1,N/4-1，则有：令 则 14042120414012214022244244NllkNkNNllkNNlklNNllkNWlxWWlxWlxWlxkA 14,1,04414NkWxkCkNNo（2.242）14,1,024414NkWxkDkNNo（2.243）14,1,02NkhDWkCkAkN（2.244）14,1,042NkkDWkCNkAkN（2.245）同样，令则有：14,1,014414NkWxkEkNNo 14,1,034414NkWxkFkNNo（2.246）（2.247）14,1,0
29、2NkkFWkEkBkN（2.248）14,1,042NkkFWkENkBkN（2.249）对于一个N=8的序列，此时的C(k)、D(k)、E(k)和F(k)均已为两点的序列，无需再分，此时有 图2.75 FFT时间抽取算法信号流图（N=8）C(0)=x(0)+x(4)，E(0)=x(1)+x(5)C(1)=x(0)-x(4)，E(1)=x(1)-x(5)D(0)=x(2)+x(6)，F(0)=x(3)+x(7)D(1)=x(2)-x(6)，F(1)=x(3)-x(7)在FFT的整个运算过程中，每两个等式的运算过程可以用一个形似蝴蝶结的“X”形结构图来表示，八个等式对应于四个蝶形结构，因此这种
30、信号流程图称为FFT的蝶形运算流程图，将这种运算的基本单元称为蝶形运算单元(butterfly computation)。图2.76 蝶形运算单元 时间抽取算法的规律：1.分级运算：将N个点的序列逐次对分，直至分到N/2个两个点的序列为止。上下节点p,q间的距离为：p-q=2m图2.77 8点FFT时间抽取算法信号流图 2.蝶形运算单元组 每一级上的N/2个蝶形单元可分为若干组，称之为蝶形运算单元组，每一组中的蝶形单元有着相同的结构和Wr因子分布。每级的蝶形单元组数目不同，第m级的组数为N/2m+1，m=0,1,m+13.Wr因子的分布 Wr因子分布的一般规律为：其中m为级次。12,1,0,12mrrWm4.数据排列顺序(data ordering)从图2.77可见，变换后的输出序列X(k)按正序排列，但在输入端序列的排列次序不是原来的自然顺序，而变成了0，4，2，6，1，5，3，7。图2.78 数据整序方法（a）奇偶分解整序 （b）码位倒置整序
