数据采集与处理

1、数据采集与处理技术电子科技大学空天科学技术研究院数据采集与应用数据采集与应用数据处理部分数据处理部分数据采集与处理技术电子科技大学空天科学技术研究院 还原采集到的电信号的物理意义还原采集到的电信号的物理意义 消除数据中的干扰信号消除数据中的干扰信号 分析计算数据的内在特征分析计算数据的内在特征把把A/D转换的数字量变换为带有工程单位具有实际物理意义的数转换的数字量变换为带有工程单位具有实际物理意义的数字量，便于进一步的处理和使用字量，便于进一步的处理和使用采用各种方法（如去除趋势项、滤波等）最大限度地消除混入数采用各种方法（如去除趋势项、滤波等）最大限度地消除混入数据中的干扰信号，以保证采集数

2、据的精度据中的干扰信号，以保证采集数据的精度对数据进行变换加工（如傅立叶变换）或相互运算（如相关函数）对数据进行变换加工（如傅立叶变换）或相互运算（如相关函数）得到能表达数据内在特征的二次数据得到能表达数据内在特征的二次数据数据采集与处理技术电子科技大学空天科学技术研究院 按处理方式划分按处理方式划分 按处理性质划分按处理性质划分在线处理：采集数据的同时对有限的数据做简单、基本的处理在线处理：采集数据的同时对有限的数据做简单、基本的处理离线处理：不受处理时间的限制，可对采集后的数据做复杂处理离线处理：不受处理时间的限制，可对采集后的数据做复杂处理预处理：如标度变换、去除趋势项等预处理：如标度变

3、换、去除趋势项等二次处理：对采集后的数据进行各种数学运算如傅立叶变换、相二次处理：对采集后的数据进行各种数学运算如傅立叶变换、相关运算等关运算等数据采集与处理技术电子科技大学空天科学技术研究院 数据采集与处理技术（第二版），马明建，西安交通大学出版社，数据采集与处理技术（第二版），马明建，西安交通大学出版社，2005.9 实用实用FFT信号处理技术，侯朝焕、阎世尊、蒋银林，海洋出版社，信号处理技术，侯朝焕、阎世尊、蒋银林，海洋出版社，1990.10 数字信号处理理论、算法与实现（第二版），胡广书，清华大学数字信号处理理论、算法与实现（第二版），胡广书，清华大学出版社，出版社，2003.8 虚拟

4、仪器及其在机械测试中的应用，柏林、王见、秦树人，科学虚拟仪器及其在机械测试中的应用，柏林、王见、秦树人，科学出版社，出版社，2007.9数据采集与处理技术电子科技大学空天科学技术研究院MATLAB (Matrix Laboratory) 是是Mathworks公司于上世纪公司于上世纪80年代年代推出的一套高性能的数值计算和可视化软件。它将一个优秀软件的推出的一套高性能的数值计算和可视化软件。它将一个优秀软件的易用性与可靠性、通用性与专业性易用性与可靠性、通用性与专业性 、一般目的应用与高深科学技术、一般目的应用与高深科学技术应用进行了有机的结合，成功地应用于各工程学科的研究领域。应用进行了有机

5、的结合，成功地应用于各工程学科的研究领域。 语言简洁紧凑，语法要求不严格，程序设计自由度较大语言简洁紧凑，语法要求不严格，程序设计自由度较大 运算符、库函数丰富运算符、库函数丰富 图形功能强大图形功能强大 界面友好，编程效率高界面友好，编程效率高数据采集与处理技术电子科技大学空天科学技术研究院 可扩展性强，应用灵活的兼容与接口功能可扩展性强，应用灵活的兼容与接口功能 具有丰富的工具箱：主工具箱、信号处理工具箱、控制系统工具具有丰富的工具箱：主工具箱、信号处理工具箱、控制系统工具箱、图像处理工具箱、通讯工具箱、系统辨识工具箱、神经网络箱、图像处理工具箱、通讯工具箱、系统辨识工具箱、神经网络工具箱

6、、金融工具箱、符号计算工具箱等工具箱、金融工具箱、符号计算工具箱等主界面及各部分功能主界面及各部分功能常用标点功能常用标点功能 “%” 注释当前行注释当前行 “” 当指令较长时，进行续行当指令较长时，进行续行数据采集与处理技术电子科技大学空天科学技术研究院路径路径窗口窗口历史历史记录记录程序程序窗口窗口指令指令窗口窗口数据采集与处理技术电子科技大学空天科学技术研究院 “,” 用于显示计算结果之间的指令分隔；数组元素间分隔用于显示计算结果之间的指令分隔；数组元素间分隔 “;” 用于不显示计算结果指令的用于不显示计算结果指令的“结尾结尾”标志；用于不显示计标志；用于不显示计算结果之间的指令分隔；数

7、组行间分隔算结果之间的指令分隔；数组行间分隔 “:” 生成一维数组；用于单下标索引时表示将数组中全部元素生成一维数组；用于单下标索引时表示将数组中全部元素构成一列；用于多下标索引时表示取数组中该维上全部数据构成一列；用于多下标索引时表示取数组中该维上全部数据常用控制指令功能常用控制指令功能 “clc” 清除指令窗中显示内容清除指令窗中显示内容 “clf” 清除图形窗清除图形窗 “clear” 清除工作空间中保存的变量清除工作空间中保存的变量数据采集与处理技术电子科技大学空天科学技术研究院 “edit” 打开打开m文件编辑器或指定函数的源文件文件编辑器或指定函数的源文件 “exit / quit

8、” 退出退出MATLAB内存变量查询功能内存变量查询功能命令：命令：who 或或 whos作用都是列出在作用都是列出在matlab工作空间中已经驻留的变量名清单，不同的工作空间中已经驻留的变量名清单，不同的是是whos在给出驻留变量的同时，还给出他们的维数性质在给出驻留变量的同时，还给出他们的维数性质 whoYour variables are:A B whos Name Size Bytes Class Attributes A 2x2 32 double B 1x6 12 char 数据采集与处理技术电子科技大学空天科学技术研究院文件操作函数文件操作函数 保存数据文件保存数据文件 save

9、函数函数 save 文件名文件名 save 文件名文件名 变量名变量名 save 文件名文件名 变量名变量名 -ascii 读取数据文件读取数据文件 load函数函数 load 文件名文件名 load 文件名文件名 变量名（限于变量名（限于“.mat”文件）文件）注：注： 缺省为缺省为“.mat”文件，文件， “ -ascii” 为文本文件为文本文件数据采集与处理技术电子科技大学空天科学技术研究院 二维图形绘制核心函数：二维图形绘制核心函数：plot绘图功能绘图功能 plot( x ,s) plot( x, y, s) plot(x1,y1,s1, x2,y2,s2, )s是用来指定线形、是用

10、来指定线形、色彩及数据点形的色彩及数据点形的选项字符串选项字符串 线形：线形：实线（实线（-）、）、 虚线（虚线（:）、点划线（）、点划线（-.）、双划线（）、双划线（-） 色彩：色彩：蓝（蓝（b）、）、 绿（绿（g）、）、 红（红（r）、）、 青（青（c） 品红（品红（m）、黄（）、黄（y）、）、 黑（黑（k）、）、 白（白（w）数据采集与处理技术电子科技大学空天科学技术研究院 各种标注函数：各种标注函数： title(图形标题图形标题) xlabel( x轴坐标名称轴坐标名称) ylabel( y轴坐标名称轴坐标名称) text(x,y,文字说明文字说明) 点形：点形：实心黑点（实心黑点（

11、.）、）、 十字符（十字符（+）、）、 八线符（八线符（*） 朝上三角（朝上三角（）、朝下三角（）、朝下三角（v）、朝左三角符（）、朝左三角符（）、菱形（）、菱形（d）、）、 六角星符（六角星符（h） 空心圆圈（空心圆圈（o）、五角星符（）、五角星符（p）、方块（）、方块（s） 叉字（叉字（x）数据采集与处理技术电子科技大学空天科学技术研究院 坐标轴设定函数：坐标轴设定函数：axis axis(xmin,xmax,ymin,ymax)设定坐标轴最大和最小值设定坐标轴最大和最小值 axis(auto)坐标轴缺省状态坐标轴缺省状态 axis(square)正方形坐标系正方形坐标系 axis(equ

12、al)纵横轴采用等长刻度纵横轴采用等长刻度 axis(on)显示坐标系统显示坐标系统 axis(off)关闭坐标系统关闭坐标系统 图例函数：图例函数：legend legend(图例图例1说明说明,图例图例2说明说明, )数据采集与处理技术电子科技大学空天科学技术研究院 其它常用函数其它常用函数 grid on网格开启网格开启 grid off网格关闭网格关闭 hold on在当前图形上继续画曲线使能在当前图形上继续画曲线使能 hold off在当前图形上继续画曲线关闭在当前图形上继续画曲线关闭 figure新建图形新建图形 绘图区分割函数：绘图区分割函数： subplot subplot(m

13、,n,p) 该函数将当前图形窗口分割成该函数将当前图形窗口分割成mn个绘图区域，即每行个绘图区域，即每行n个，共个，共m行区域，且选定第行区域，且选定第p个区为当前区。个区为当前区。数据采集与处理技术电子科技大学空天科学技术研究院数据产生函数数据产生函数程序文件（程序文件（m文件）文件） “正弦信号正弦信号” sin函数函数 “方波信号方波信号” square函数函数 “调频余弦信号调频余弦信号” chirp函数函数 “sinc信号信号” sinc函数函数 “随机噪声随机噪声” randn函数函数一种要实现的命令写在一个以一种要实现的命令写在一个以m为文件扩展名的文件中，然后由为文件扩展名的文

14、件中，然后由MATLAB系统解释执行系统解释执行 数据采集与处理技术电子科技大学空天科学技术研究院 脚本文件：脚本文件： 函数文件函数文件 一种纯粹的命令指令的集合体，它没有输入输出参数；脚本一种纯粹的命令指令的集合体，它没有输入输出参数；脚本M文文件中的变量为全局变量，在执行过程中，脚本文件可以使用当前件中的变量为全局变量，在执行过程中，脚本文件可以使用当前工作区中的变量，而脚本文件中产生的变量，也保留在工作区内工作区中的变量，而脚本文件中产生的变量，也保留在工作区内 一般作为子程序来用，经常被脚本一般作为子程序来用，经常被脚本M文件调用，甚至也可以被其文件调用，甚至也可以被其它函数它函数M

15、文件调用，其格式为：文件调用，其格式为：function z = add (x,y) 文件必须以文件必须以function开头，函数名必须与将来保存的开头，函数名必须与将来保存的m文件名一文件名一致，可有输入输出参数，且变量皆为局部变量致，可有输入输出参数，且变量皆为局部变量数据采集与处理技术电子科技大学空天科学技术研究院画出衰减振荡曲线画出衰减振荡曲线及其它的包络及其它的包络 t取值范围取值范围 。例例1/3sin3tyet/30tye0,4024681012-1-0.8-0.6-0.4-0.200.20.40.60.81tamplitudey=exp(-t/3)*sin(3t) y y0-

16、y0数据采集与处理技术电子科技大学空天科学技术研究院程序清单：程序清单：数据采集与处理技术电子科技大学空天科学技术研究院例例2 在某工程实际在某工程实际应用中，有一信号的应用中，有一信号的主要频率成分是主要频率成分是50Hz和和300Hz的正弦的正弦信号，该信号被一白信号，该信号被一白噪声污染，现对该信噪声污染，现对该信号进行采样，采样频号进行采样，采样频率为率为1000Hz，通过，通过傅立叶变换对其频率傅立叶变换对其频率成份进行分析。成份进行分析。数据采集与处理技术电子科技大学空天科学技术研究院程序清单：程序清单：数据采集与处理技术电子科技大学空天科学技术研究院把把A/D转换的数字量变换为带

17、有工程单位的数字量的变换称为标度转换的数字量变换为带有工程单位的数字量的变换称为标度变换（工程变换）。标度变换有多种形式，取决于被测物理量所用变换（工程变换）。标度变换有多种形式，取决于被测物理量所用的传感器或变送器的类型。的传感器或变送器的类型。线性参数的标度变换线性参数的标度变换当被测物理量与传感器或仪表的输出之间呈线性关系时采用线性变当被测物理量与传感器或仪表的输出之间呈线性关系时采用线性变换，其公式为：换，其公式为：m000m0Y -YY=Y +(X-N )N-N数据采集与处理技术电子科技大学空天科学技术研究院Y 标度变换后所得到的被测量的实际值标度变换后所得到的被测量的实际值X 被测

18、量实际值被测量实际值Y所对应的所对应的A/D转换后的数字量转换后的数字量Y0 被测量量程的下限被测量量程的下限Ym 被测量量程的上限被测量量程的上限N0 Y0对应的对应的A/D转换后的数字量转换后的数字量Nm Ym对应的对应的A/D转换后的数字量转换后的数字量适用于被测物理量与传感器或变送器的输出之间的关系为非线性但适用于被测物理量与传感器或变送器的输出之间的关系为非线性但可用解析式表达的情况。可用解析式表达的情况。其中：其中：非线性参数的标度变换非线性参数的标度变换 公式变换法公式变换法数据采集与处理技术电子科技大学空天科学技术研究院如：如：在流量测量中，从差压变送器来的信号在流量测量中，从

19、差压变送器来的信号P与实际流量与实际流量Q成平方成平方关系，即关系，即 。式中。式中K为刻度系数，与流体的性质和节流装为刻度系数，与流体的性质和节流装置尺寸有关。则标度变换公式为置尺寸有关。则标度变换公式为m000m0Q-QY = Q+(X -N)N-NQ=KPY 被测流量经标度变换的实际值被测流量经标度变换的实际值X 被测流量实际值被测流量实际值Y所对应的所对应的A/D转换后的数字量转换后的数字量Q0 被测量量程的下限被测量量程的下限Qm 被测量量程的上限被测量量程的上限N0 Q0对应的对应的A/D转换后的数字量转换后的数字量Nm Qm对应的对应的A/D转换后的数字量转换后的数字量其中：其中

20、：数据采集与处理技术电子科技大学空天科学技术研究院适用于被测物理量与传感器或变送器的输出之间的关系为非线性且适用于被测物理量与传感器或变送器的输出之间的关系为非线性且无法用解析式表达的情况。寻找多项式的方法有最小二乘法、代数无法用解析式表达的情况。寻找多项式的方法有最小二乘法、代数插值法等。插值法等。 多项式变换法多项式变换法代数插值法原理代数插值法原理已知被测量已知被测量 与传感器的输出值与传感器的输出值 在在 个相异点个相异点处的函数值为处的函数值为0011(),(),()nnf xyf xyf xy012naxxxxb( )yf xx1n数据采集与处理技术电子科技大学空天科学技术研究院可

21、解得可解得 。用一个次数不超过用一个次数不超过 的代数多项式的代数多项式逼近函数逼近函数 ，使，使 在点在点 处满足式：处满足式：( )( )niiiP xf xy1110( )nnnnnP xa xaxa xa( )yf xnixnP该式中待定系数该式中待定系数 ，共有，共有 个，因此只要用个，因此只要用 个已个已知的知的 和和 代入该式得方程组代入该式得方程组:01,na aa1n1nix(0,1,2, )iy in( )nP x10101000111111011110nnnnnnnnnnnnnnnna xaxa xaya xaxa xaya xaxa xay数据采集与处理技术电子科技大学

22、空天科学技术研究院假设采样数据为假设采样数据为 ，并设采样时间间隔为，并设采样时间间隔为t=1对于实际采集到的数据，由于各种因素的干扰影响，往往会偏离基对于实际采集到的数据，由于各种因素的干扰影响，往往会偏离基线，甚至偏离基线的大小还会随时间变化。偏离基线随时间变化的线，甚至偏离基线的大小还会随时间变化。偏离基线随时间变化的整个过程称为信号的趋势项。常用最小二乘法消除趋势项。整个过程称为信号的趋势项。常用最小二乘法消除趋势项。最小二乘法消除趋势项的原理最小二乘法消除趋势项的原理 (1,2,3, )kxkn设一个多项式函数为：设一个多项式函数为：2201 (1,2,3, )mkmxaaka ka

23、 kkn数据采集与处理技术电子科技大学空天科学技术研究院满足满足E有极值的条件为：有极值的条件为：确定函数确定函数 的各待定系数的各待定系数 ，使得函数，使得函数 与离散与离散数据数据 的误差平方和最小，即的误差平方和最小，即依次取依次取E对求偏导，可以产生一个对求偏导，可以产生一个m+1元线性方程组：元线性方程组：kx (0,1, )jajmkxkx22110()()knnmjkjkkkjExxa kx 102()0 (0,1,)nmijjkkjiEka kxima1010 (0,1,)nmnj iijkkjka kx kim数据采集与处理技术电子科技大学空天科学技术研究院解方程组，求出解方

24、程组，求出m+1个待定系数个待定系数 。上面各式中，。上面各式中，m为设定的多项式阶次。为设定的多项式阶次。解得：解得： (0,1, )jajm000110 nnkkka kx k讨论：讨论： 当当m=0时为常数趋势项，有时为常数趋势项，有011 nkkaxn则消除常数趋势项的公式为：则消除常数趋势项的公式为：0 (1,2,3, )kkkkyxxxakn数据采集与处理技术电子科技大学空天科学技术研究院此时可求解方程得此时可求解方程得 的值的值00011112011110 0nnnkkkknnnkkkka ka kx ka ka kx k 当当m=1时为线性趋势项，有时为线性趋势项，有01aa和

25、则消除线性趋势项的公式为：则消除线性趋势项的公式为：01(+) (1,2,3, )kkkkyxxxaa kkn数据采集与处理技术电子科技大学空天科学技术研究院说明：说明： 当当m2时为曲线趋势项，用同样的方法也可以求出时为曲线趋势项，用同样的方法也可以求出 的各个系的各个系数。代入数。代入 即可消除趋势项。即可消除趋势项。 多项式阶次的选择要根据物理过程的情况而定。计算阶次高，当多项式阶次的选择要根据物理过程的情况而定。计算阶次高，当然拟合的精度高，即所求的趋势项精度也高，但计算量相应也复然拟合的精度高，即所求的趋势项精度也高，但计算量相应也复杂。在实际的数据处理中，通常取杂。在实际的数据处理

26、中，通常取m=13。 (1,2,3, )kkkyxxknkx 消除趋势项的工作要特别谨慎，只有当采样数据中有明显的趋势消除趋势项的工作要特别谨慎，只有当采样数据中有明显的趋势项时才需考虑消除。项时才需考虑消除。数据采集与处理技术电子科技大学空天科学技术研究院举例：举例：仿真一含有随机噪声的正弦信号，其频率为仿真一含有随机噪声的正弦信号，其频率为7Hz，加入一线，加入一线性趋势项性趋势项 ，试用最小二乘法去除该趋势项。，试用最小二乘法去除该趋势项。MATLAB中关键实现函数：中关键实现函数： 最小二乘法定阶数的多项式拟合：最小二乘法定阶数的多项式拟合：a= polyfit ( t,x,m) 。用

27、最小二。用最小二乘法对输入数据乘法对输入数据 t 和和 x 用用 m 阶多项式进行逼近，函数返回一个阶多项式进行逼近，函数返回一个长度为长度为m+1的多项式的系数，并存到向量中。的多项式的系数，并存到向量中。求多项式的值：求多项式的值：b=polyval (a,t)。表示以输入变量。表示以输入变量 t 向量，计算向量，计算系数向量系数向量 a 的多项式的值，计算结果存到向量的多项式的值，计算结果存到向量 b 中。中。( )1 5 x tt 数据采集与处理技术电子科技大学空天科学技术研究院程序清单：程序清单：数据采集与处理技术电子科技大学空天科学技术研究院00.10.20.30.40.50.60

28、.70.80.9102468tSignal With Trend sin+trendtrend00.10.20.30.40.50.60.70.80.91-2-1012tSignal Without Trend 运行结果：运行结果：数据采集与处理技术电子科技大学空天科学技术研究院波形分析包括幅值参数和一些由幅值参数演化而来的无量纲指标。波形分析包括幅值参数和一些由幅值参数演化而来的无量纲指标。均值和绝对均值均值和绝对均值均值指信道中的直流分量，是信号在分析区间内的算术平均。均值指信道中的直流分量，是信号在分析区间内的算术平均。01( )Txxtd tT11()nkxxkN绝对均值指信号绝对值的算

29、术平均，是信号强度的平均。绝对均值指信号绝对值的算术平均，是信号强度的平均。01( )Txxtd tT11()nkxxkN数据采集与处理技术电子科技大学空天科学技术研究院均方值（平均功率）和有效值（均方根值）均方值（平均功率）和有效值（均方根值）均方值是指信号幅值在分析区间内的平方均值。均方值是指信号幅值在分析区间内的平方均值。201( )Trxx t dtT211()nrkxxkN有效值是指均方值的平方根，具有幅值量纲，是反映确定性信号作有效值是指均方值的平方根，具有幅值量纲，是反映确定性信号作用强度的主要时域参数。用强度的主要时域参数。201( )TRMSxxt dtT211()nRM S

30、kxxkN数据采集与处理技术电子科技大学空天科学技术研究院峰值（单峰值）和双峰值峰值（单峰值）和双峰值峰值指分析区间内出现的最大幅值，它反映了信号的瞬时最大作用峰值指分析区间内出现的最大幅值，它反映了信号的瞬时最大作用强度。强度。max( )min( )ppxx kx k双峰值又称峰双峰值又称峰-峰值（峰值（peak-peak），是指正负峰值（最小值）之间），是指正负峰值（最小值）之间的差，它不仅反映了信号的瞬时作用强度，还反映了信号幅值的变的差，它不仅反映了信号的瞬时作用强度，还反映了信号幅值的变化范围和偏离中心位置的情况。化范围和偏离中心位置的情况。max( )pxx k数据采集与处理技术

31、电子科技大学空天科学技术研究院无量纲指标无量纲指标无量纲指标可以从不同侧面反映信号波形形状及特性。主要包括：无量纲指标可以从不同侧面反映信号波形形状及特性。主要包括：波形因素波形因素K (kurtosis factor)RMSxKx峰值因素峰值因素C (crest factor)pRMSxCx脉冲因素脉冲因素I (impulse factor)pxIx裕度因素裕度因素L (clearance factor)prxLx数据采集与处理技术电子科技大学空天科学技术研究院相关分析在检测系统、控制系统、通信系统等领域广为应用，它主相关分析在检测系统、控制系统、通信系统等领域广为应用，它主要解决信号内部联

32、系问题，信号与信号之间的相似性问题。要解决信号内部联系问题，信号与信号之间的相似性问题。定义定义对于两信号对于两信号 ，相关函数定义为：，相关函数定义为：( )( )x ny n、对于实际应用的物理信号是因果性的，即当对于实际应用的物理信号是因果性的，即当 时，时， 恒为恒为0则则1()lim( )()21NxyNnNrmx n y nmN 0n ( )x n101()lim( )()NxyNnrmx n y nmN 数据采集与处理技术电子科技大学空天科学技术研究院性质性质相关函数描述了两个信号在延时相关函数描述了两个信号在延时m时的相似程度。时的相似程度。 自相关函数自相关函数 的性质：的性

33、质：( )xr m101()lim( ) ()NxNnrmx n x nmN 当当 时，称为自相关函数：时，称为自相关函数：( )( )x ny n= 为偶函数且当为偶函数且当 时取得最大值；时取得最大值； 若若 是周期为是周期为T的周期信号，则自相关函数也为同样周期的的周期信号，则自相关函数也为同样周期的周期函数；周期函数；0m ( )x n数据采集与处理技术电子科技大学空天科学技术研究院 若若 不包含任何周期分量，则当不包含任何周期分量，则当m越大时信号自相关性越低，越大时信号自相关性越低，m趋于无穷时，趋于无穷时， 与与 彼此不相关。彼此不相关。 互相关函数互相关函数 的性质：的性质：(

34、)xyrm 既不是奇函数也不是偶函数，；既不是奇函数也不是偶函数，； 若若 具有相同频率的周期成分，则具有相同频率的周期成分，则 中包含有该中包含有该频率的周期成分，在频率的周期成分，在 时，该频率的周期仍然存在；时，该频率的周期仍然存在； 在在m = 0时，时， 一般不会得到最大值。一般不会得到最大值。( )()xyyxrmrm( )x n( )( )x ny n、()xyrmm ()xyrm( )x n()x nm数据采集与处理技术电子科技大学空天科学技术研究院相关函数的应用相关函数的应用 检查混淆在随机信号中的确定性周期信号检查混淆在随机信号中的确定性周期信号例：例：仿真一混有频率为仿真

35、一混有频率为5Hz正弦信号的随机噪声数据，设采样频率正弦信号的随机噪声数据，设采样频率为为1000Hz，采样时间为，采样时间为5秒，试用相关函数法检测出该正弦信号。秒，试用相关函数法检测出该正弦信号。数据采集与处理技术电子科技大学空天科学技术研究院00.511.522.533.544.55-4-2024tOrigin Signal-5-4-3-2-1012345-0.500.511.5tCorrelation Function运行结果：运行结果：数据采集与处理技术电子科技大学空天科学技术研究院 测定深埋地下的输油管道裂损的位置测定深埋地下的输油管道裂损的位置示意图如下：示意图如下：数据采集与处

36、理技术电子科技大学空天科学技术研究院漏损处视为传播声源，两侧管道分别放置传感器，因为放置传感器漏损处视为传播声源，两侧管道分别放置传感器，因为放置传感器的两点距漏损处的距离不相等，放漏油的音响传至两个传感器时就的两点距漏损处的距离不相等，放漏油的音响传至两个传感器时就有时差。在互相关图上有时差。在互相关图上 处有最大值，这个处有最大值，这个 就是时差。由此就是时差。由此即可确定漏损处的位置：即可确定漏损处的位置：dd12ddV式中：式中： 两传感器的中点距离漏油处的距离；两传感器的中点距离漏油处的距离； 声响通过管道的传播速度。声响通过管道的传播速度。dV数据采集与处理技术电子科技大学空天科学

37、技术研究院相关函数的估计相关函数的估计下面以自相关函数为例说明相关函数的估计问题。下面以自相关函数为例说明相关函数的估计问题。 直接估计法直接估计法若若 的总点数的总点数N为有限值，则求为有限值，则求 估计值的一种方法是：估计值的一种方法是：( )x n101( )( )()NxNNnrmxn xnmN( )xr m由于由于 只有只有N个值，因此对于每个固定的延迟个值，因此对于每个固定的延迟m，可以利用的数，可以利用的数据只有据只有 个，所以实际计算时上式变为：个，所以实际计算时上式变为：( )x nNm101( )( )()NmxNNnrmxn xnmN 数据采集与处理技术电子科技大学空天科

38、学技术研究院 对对 作逆变换，得作逆变换，得上式不是上式不是 的无偏估计，修改为无偏估计如下：的无偏估计，修改为无偏估计如下： 间接估计法（基于间接估计法（基于FFT的快速计算）的快速计算）( )xr m101( )( )()NmxNNnrmxn xnmNm 对对 补补 个个0，得，得 ，对，对 做做DFT得得 ，其，其中中 求求 的幅平方，然后除以的幅平方，然后除以 得得( )NxnN2( )Nxn2( )Nxn2( )NXk0,1,21kN2( )NXkN221( )NXkN221( )NXkN0( )rm数据采集与处理技术电子科技大学空天科学技术研究院 并不简单地等于并不简单地等于 ，而

39、是等于将，而是等于将 中中的部分向右平移的部分向右平移 点后形成的新序列，如下图所示。因此，只需点后形成的新序列，如下图所示。因此，只需将所求的将所求的 稍作修改即可得到稍作修改即可得到 。0( )rm( )r m( )r m(1)0Nm2N0( )rm( )r m证明详见：证明详见：数字信号处理理论算法与实现数字信号处理理论算法与实现,胡广书胡广书,清华大学出版社清华大学出版社数据采集与处理技术电子科技大学空天科学技术研究院补充内容：基于补充内容：基于FFT的互相关函数的快速估计法的互相关函数的快速估计法设设 的点数都为有限值的点数都为有限值N ，则，则 的估计值是：的估计值是：( )( )

40、x ny n、101()( )()NmxyNNnrmxn ynmN ( )xyrm有偏估计：有偏估计：101( )( )()NmxyNNnrmxn ynmNm 无偏估计：无偏估计：基于基于FFT的互相关函数的快速估计算法步骤如下：的互相关函数的快速估计算法步骤如下：数据采集与处理技术电子科技大学空天科学技术研究院 将将 后面的后面的 个数据移到最前面得调整位置后的个数据移到最前面得调整位置后的 对对 作逆变换，得作逆变换，得 对对 分别补分别补 个个0，得，得 ，对其做，对其做DFT得得 ，其中，其中 求求 ( )( )Nxny n、N22( )( )NNxnyn、22( )( )NNXkYk

41、、0,1,21kN22( )( )NNXk Yk1N ( )xyrm22( )( )NNXk Yk( )xyrm( )xyrm 对对 进行系数修正，有偏估计为进行系数修正，有偏估计为 ，无偏估计为，无偏估计为( )xyrm1( )xyrmN1( )xyrmNm数据采集与处理技术电子科技大学空天科学技术研究院程序清单：程序清单：数据采集与处理技术电子科技大学空天科学技术研究院周期连续信号周期连续信号傅立叶级数傅立叶级数若周期函数若周期函数 的周期为的周期为T且满足狄里赫利条件，则此周期函数可且满足狄里赫利条件，则此周期函数可表示为傅立叶级数的三角函数形式：表示为傅立叶级数的三角函数形式：2021

42、( )TTax t dtT式中：式中：00001( )(cossin) (2)nnnx taantbntT其中2022( ) c o sTTnax tntd tT2022( ) sinTTnbx tntd tT( )x t数据采集与处理技术电子科技大学空天科学技术研究院进一步化简：进一步化简：式中：式中：001( )cos() nnnx taAnt22nnnAabnnnba rctga 也可表示为复指数形式（由欧拉公式即可得到）：也可表示为复指数形式（由欧拉公式即可得到）：0( ) jntnnx tC e式中：式中：0221( )nTjntjTnnCx t edtCeT12nnCAnnnba

43、rctga 数据采集与处理技术电子科技大学空天科学技术研究院对于如图的周期方波函数，一个周期内的数学表达为：对于如图的周期方波函数，一个周期内的数学表达为：0004A11( )(sinsin 3sin 5)35x ttttA 02 ( )-A -20tTx tTt 因此可表示为：因此可表示为：0002(1cos)nnaaAbnn代入傅立叶系数公式得：代入傅立叶系数公式得：数据采集与处理技术电子科技大学空天科学技术研究院数据采集与处理技术电子科技大学空天科学技术研究院周期信号的特点：周期信号的特点： 周期信号的频谱是由无限多条离散谱线组成，每一条谱线（单周期信号的频谱是由无限多条离散谱线组成，每

44、一条谱线（单边谱）代表一个谐波分量；边谱）代表一个谐波分量； 各次谐波的频率只能是基波频率的整数倍；各次谐波的频率只能是基波频率的整数倍； 谱线的高度表示了相应谐波分量的幅值大小。对于工程常见的谱线的高度表示了相应谐波分量的幅值大小。对于工程常见的周期信号，其谐波幅值的总趋势是随着谐波次数的增高而减小。周期信号，其谐波幅值的总趋势是随着谐波次数的增高而减小。当谐波次数无限增高时，其幅值趋于当谐波次数无限增高时，其幅值趋于0.由图可知：方波的幅频谱只包含基波和奇次谐波的频率分量，且谐由图可知：方波的幅频谱只包含基波和奇次谐波的频率分量，且谐波的幅值以波的幅值以 的规律收敛；相频谱各次谐波的初相位

45、皆为的规律收敛；相频谱各次谐波的初相位皆为0。1 n数据采集与处理技术电子科技大学空天科学技术研究院非周期连续信号非周期连续信号傅立叶变换傅立叶变换000221( ) =( ) ( )21 ( )2TjntjntjntTnnnjtjtjtjtx tC ex t edteTdx t edtex t edted 思路：非周期信号可以看成是周期思路：非周期信号可以看成是周期T趋于无穷大的周期信号。当周趋于无穷大的周期信号。当周期期T延拓时，区间从延拓时，区间从 趋于趋于 ，频谱的频率间隔变化方，频谱的频率间隔变化方式式 ，离散的，离散的 变为连续的变为连续的 ，展开式的叠，展开式的叠加关系变为积分关

46、系。此时周期信号的复指数展开形式变为：加关系变为积分关系。此时周期信号的复指数展开形式变为：,22TT, 02Td0n数据采集与处理技术电子科技大学空天科学技术研究院记：记：( )j tXx t edt21( )()()2jtjftx tXedXfedf称称 为为 的傅立叶变换，的傅立叶变换， 为为 的逆傅立叶变换，两的逆傅立叶变换，两者互称为傅立叶变换对。者互称为傅立叶变换对。 X( )x t( )x t X例：例：求矩形窗函数的频谱，矩形窗定义如下：求矩形窗函数的频谱，矩形窗定义如下：1 2( )0 2tTw ttT2( )jftXfx t edt或或数据采集与处理技术电子科技大学空天科学

47、技术研究院解：解：2222( ) 1 2sin sinc()jftTjftTjfTjfTWfw t edtedteejffTTfTTfT其中，其中，sinc(x)函数是偶函数以函数是偶函数以 为周期并随为周期并随x的增加而做衰减振荡，的增加而做衰减振荡，在在 处其值为处其值为0。2 (1, 2,)nn 数据采集与处理技术电子科技大学空天科学技术研究院0210()( )jf tFXffx t e式中式中 为常数，表明：若频谱沿频率轴平移一个常值为常数，表明：若频谱沿频率轴平移一个常值 ，对应的，对应的时域函数将乘时域函数将乘0f傅立叶变换两条重要性质：傅立叶变换两条重要性质： 频移性质频移性质0

48、f02ff te 卷积性质卷积性质( )( )()()( )( )()()F x ty tXf YfF x t y tXfYf式中：式中：( )( )( )()x ty txy td该性质表明：时域乘积对应频域卷积，时域卷积对应频域乘积。该性质表明：时域乘积对应频域卷积，时域卷积对应频域乘积。数据采集与处理技术电子科技大学空天科学技术研究院离散数字信号离散数字信号离散傅立叶变换（离散傅立叶变换（DFT）正变换和逆变换的公式：正变换和逆变换的公式：21100( ) ( )( )( )NNjnknkNNnnX kDFT x nx n ex n W2110011( )( )( )( )NNjnknk

49、NNkkx nIDFT X kX k eX k WNN在计算机上用在计算机上用FFT（快速傅立叶变换）算法计算离散傅立叶变换的（快速傅立叶变换）算法计算离散傅立叶变换的结果。其复数的乘法运算量由结果。其复数的乘法运算量由 提高到提高到2N2log2NN数据采集与处理技术电子科技大学空天科学技术研究院信号的截断及能量泄露效应信号的截断及能量泄露效应傅立叶变换是研究整个时间傅立叶变换是研究整个时间域与频率域的关系。然而计域与频率域的关系。然而计算机处理实际数据时，只能算机处理实际数据时，只能处理有限个数据，相当于将处理有限个数据，相当于将无限长的信号进行了截断，无限长的信号进行了截断，即将无限长的

50、信号乘以一个即将无限长的信号乘以一个窗函数。下面以余弦信号为窗函数。下面以余弦信号为例说明：例说明：数据采集与处理技术电子科技大学空天科学技术研究院由频域卷积定理知截断后的余弦信号的频谱应为原始无限长余弦信由频域卷积定理知截断后的余弦信号的频谱应为原始无限长余弦信号的频谱与矩形窗函数频谱的卷积。其频谱发生了畸变，原来集中号的频谱与矩形窗函数频谱的卷积。其频谱发生了畸变，原来集中在在 处的能量被分散到两个较宽的频带中，称为能量泄露。处的能量被分散到两个较宽的频带中，称为能量泄露。0f窗函数窗函数时间域中可采用不同的窗函数减小泄露的影响。窗函数选取原则：时间域中可采用不同的窗函数减小泄露的影响。窗

51、函数选取原则： 主瓣尽量窄以提高频率分辨率；主瓣尽量窄以提高频率分辨率； 旁瓣尽量低，以减少泄露。旁瓣尽量低，以减少泄露。两者往往不能同时满足，需根据不同的测试对象选择不同的窗函数两者往往不能同时满足，需根据不同的测试对象选择不同的窗函数数据采集与处理技术电子科技大学空天科学技术研究院从减少泄露（旁瓣要窄）从减少泄露（旁瓣要窄）的观点看，汉宁窗优于矩的观点看，汉宁窗优于矩形窗，但汉宁窗主瓣较宽，形窗，但汉宁窗主瓣较宽，频率分辨率下降。频率分辨率下降。常用的窗函数有：矩形窗（常用的窗函数有：矩形窗（boxcar）、汉宁窗（）、汉宁窗（hanning）、海明）、海明窗（窗（hamming）、三角窗

52、（）、三角窗（triang）、布莱克曼窗（）、布莱克曼窗（blackman）等。下图是矩形窗和汉宁窗的频谱对比图。等。下图是矩形窗和汉宁窗的频谱对比图。数据采集与处理技术电子科技大学空天科学技术研究院采样频率采样频率Fs、数据点数、数据点数N、采样时间、采样时间T、分析最大频率（奈奎斯特、分析最大频率（奈奎斯特频率频率fmax）频率分辨率）频率分辨率f、采样间隔、采样间隔t之间关系：之间关系：max211FsfTNtFsNftFsfT max2NffFsfN 数据采集与处理技术电子科技大学空天科学技术研究院下面以产生长度均为下面以产生长度均为n的窗函数为例说明的窗函数为例说明( )1 (0,1

53、,1)iin补充内容：窗函数公式补充内容：窗函数公式 矩形窗（矩形窗（boxcar（n） 汉宁窗（汉宁窗（hanning（n）2( )0.5 1cos() (0,1,1)iiinn数据采集与处理技术电子科技大学空天科学技术研究院2( )1 (0,1,1)iniinn 三角窗（三角窗（boxcar（n） 海明窗（海明窗（hamming（n）2( )0.540.46cos() (0,1,1)iiinn 布莱克曼窗（布莱克曼窗（blackman（n）24( )0.420.5cos() +0.08cos() (0,1,1)iiiinnn数据采集与处理技术电子科技大学空天科学技术研究院510152025

54、3035404500.20.40.60.81 Rectangularhanninghammingtriangblackman几种窗函数时域图形对比几种窗函数时域图形对比数据采集与处理技术电子科技大学空天科学技术研究院离散傅立叶变换的基本公式为：离散傅立叶变换的基本公式为：210( )( )NjnkNnX kx n e则其幅频谱为：则其幅频谱为：1( )X kfN相频谱为：相频谱为：( )X kf2( )X kfN双边谱双边谱单边谱单边谱定义定义数据采集与处理技术电子科技大学空天科学技术研究院例：例：仿真包含两个正弦（频率分别为仿真包含两个正弦（频率分别为200Hz、208Hz，幅值分别为，幅值

55、分别为1和和0.008）的信号画出其频谱，并验证加窗的效果。）的信号画出其频谱，并验证加窗的效果。数据采集与处理技术电子科技大学空天科学技术研究院运行结果：运行结果：050100150200250300350400450500-250-200-150-100-500FFT Amplitude Spectralf/Hzamp/dB Nohanning数据采集与处理技术电子科技大学空天科学技术研究院对自相关函数和互相关函数分别作傅立叶变换，得：对自相关函数和互相关函数分别作傅立叶变换，得：称称 为信号为信号 的自功率谱，的自功率谱， 为信号为信号 的互功率的互功率谱。谱。2( )( )jfmxxm

56、Pfr m e定义定义2( )( )jfmxyxymPfrm e( )xP f( )x n( )xyPf( )( )x ny n、功率谱反映了信号的功率功率谱反映了信号的功率在频率域随频率在频率域随频率f的分布，的分布，因此它们也称功率谱密度因此它们也称功率谱密度数据采集与处理技术电子科技大学空天科学技术研究院功率谱的估计功率谱的估计对于含有对于含有N个点的实际数据个点的实际数据 可以通过两种方法估计其可以通过两种方法估计其自功率谱和互功率谱：自功率谱和互功率谱：( )( )NNxnyn、 间接法：先求出相关函数，然后求其傅立叶变换得到自或互功率间接法：先求出相关函数，然后求其傅立叶变换得到自

57、或互功率谱，由于该法是通过相关函数求出，因此称为间接法。谱，由于该法是通过相关函数求出，因此称为间接法。N较小时，较小时，计算量不是太大；但当计算量不是太大；但当N较大时，在时域中计算相关函数较慢，较大时，在时域中计算相关函数较慢，是是FFT算法问世前常用的谱估计方法。算法问世前常用的谱估计方法。 直接法：直接在频域中进行估计的算法。当直接法：直接在频域中进行估计的算法。当 FFT问世后成为谱估问世后成为谱估计方法的主流。计方法的主流。数据采集与处理技术电子科技大学空天科学技术研究院211( )( )( )( )xPkX k XkX kNN 直接法（周期图法）直接法（周期图法）1( )( )(

58、 )xyPkX k YkN例：例：仿真包含两正弦（频率仿真包含两正弦（频率为为80Hz、140Hz，幅值为，幅值为1和和2）并加入随机噪声的信）并加入随机噪声的信号，用直接法估计功率谱。号，用直接法估计功率谱。00.10.20.30.40.50.60.70.80.91-6-4-20246时 间 /s幅值原 始 信 号数据采集与处理技术电子科技大学空天科学技术研究院程序清单：程序清单：数据采集与处理技术电子科技大学空天科学技术研究院数据采集与处理技术电子科技大学空天科学技术研究院结论：结论：由图可知，基本的周期图估计方法效果并不好，它的估计方由图可知，基本的周期图估计方法效果并不好，它的估计方差

59、很大，而且不满足一致性估计的条件，即方差不会随着差很大，而且不满足一致性估计的条件，即方差不会随着N的增大的增大而趋于而趋于0。改进思路：改进思路：可以将一长度为可以将一长度为N的数据的数据 分成分成L段，分别求每一段段，分别求每一段的功率谱，然后加以平均，以达到希望的目的。的功率谱，然后加以平均，以达到希望的目的。( )Nxn 直接法的改进直接法的改进1（Bartlett法，平均周期图法）法，平均周期图法）下面以自功率谱的估计为例说明，互功率谱方法类似。下面以自功率谱的估计为例说明，互功率谱方法类似。将采样数据将采样数据 分成分成L段，每段的长度都是段，每段的长度都是M，即，即N=LM，则第

60、，则第i段的功率谱为：段的功率谱为：( )Nxn数据采集与处理技术电子科技大学空天科学技术研究院把对应把对应 相加，再取平均，得到平均周期图相加，再取平均，得到平均周期图 为：为：以上例为例，将原始数据分以上例为例，将原始数据分3段，每段段，每段256个数据，得其平均周期图个数据，得其平均周期图法的功率谱过程为：法的功率谱过程为：211( )( )( )( ) (0,1)iiiixPkXk XkXkkMMM( )ixPk( )xPk21111( )( )( )LLiixxiiP kPkX kLML数据采集与处理技术电子科技大学空天科学技术研究院程序清单及运行结果：程序清单及运行结果：数据采集与