机械工程测试技术基础

1、学习要求:1数字信号处理概述2了解信号模数转换原理3掌握信号采样定理，能正确选择采样频率4了解信号截断、能量泄露、栅栏效应等现象5DFT与FFT6栅栏效应与窗函数7常用的数字信号处理算法1 数字信号处理概述 1)数字信号处理的主要研究内容 数字信号处理主要研究用数字序列来表示测试信号，并用数学公式和运算来对这些数字序列进行处理。其主要内容包括数字波形分析、幅值分析、频谱分析和数字滤波。2)测试信号数字化处理的基本步骤 物理信号对象传感器电信号信号调制电信号A/D转换数字信号计算机显示D/A转换控制3) 数字信号处理的优势 (1)用数学计算和计算机显示代替复杂的电路和机械结构(2)计算机软硬件技

2、术发展的有力推动a)多种多样的工业用计算机。 b)灵活、方便的计算机虚拟仪器开发系统2 模数(A/D)1) A/D转换过程 采样利用采样脉冲序列p(t)，从连续时间信号x(t)中抽取一系列离散样值，使之成为采样信号x(nTs)的过程 编码将离散幅值经过量化以后变为二进制数字的过程。 量化把采样信号x（nTs）经过舍入变为只有有限个有效数字的数，这一过程称为量化如4位码，只表示24=16种不同的信号幅度，这些幅度称为量化电平。当离散时间信号幅度与量化电平不相同时，就要以最接近的一个量化电平来近似它。所以经过A/D变换器后，不但时间离散化了，而且幅度也量化了，产生一个二进制流。2) A/D转换器的

3、技术指标 分辨率 用输出二进制数码的位数表示。位数越多，量化误差越 小，分辨力越高。常用有8位、10位、12位、16位等。 转换精度 转换精度LSB/2 LSB: Least Significant Bit 最低有效位转换速度 完成1次转换所用的时间,如100ms(10Hz);10us(100kHz) 模拟信号的输入范围;如5V， +/-5V，10V，+/-10V等。3 采样定理 1) 采样信号的频谱 采样过程是将采样脉冲序列p(t)与信号x(t)相乘来2) 频混现象 频域解释时域解释3)采样定理 为保证采样后信号能真实地保留原始模拟信号信息，信号采样频率必须至少为原信号中最高频率成分的2倍。

4、这是采样的基本法则，称为采样定理。fs2fmax 需要注意的是，在对信号进行采样时，满足了采样定理，只能保证不发生频率混叠，只能保证对信号的频谱作逆傅立叶变换时，可以完全变换为原时域采样信号xs(t)，而不能保证此时的采样信号能真实地反映原信号x(t)。工程实际中采样频率通常大于信号中最高频率成分的3到5倍。频混现象实验： 频混计算： FsFsFsFs频混正常Fs/2工程处理：混迭频率=Fs-信号频率A/D采样前的抗混迭滤波： 物理信号对象传感器电信号信号调理电信号A/D转换数字信号计算机显示展开低通滤波(0-Fs/2)放大4 信号的截断、能量泄漏 为便于数学处理，通常对截断的信号做周期延拓，

5、得到虚拟的无限长的信号。 当运用计算机进行测试信号处理时，不可能对无限长的信号进行测量和运算，而是取其有限的时间片段进行分析，这个信号截取过程成为信号的截断。 周期延拓后的信号与真实信号是不同的，下面我们就从数学的角度来看这种处理带来的误差情况。 设有余弦信号x(t), 用矩形窗函数w(t)与其相乘，得到截断信号: y(t) =x(t)w(t) 将截断信号谱 XT()与原始信号谱X()相比较可知，它已不是原来的两条谱线，而是两段振荡的连续谱. 原来集中在f0处的能量被分散到两个较宽的频带中去了，这种现象称之为频谱能量泄漏。能量泄漏实验： 6.5 DFT与FFT 1)离散傅立叶变换 离散傅里叶变

6、换（Discrete Fourier Transform）一词是为适应计算机作傅里叶变换运算而引出的一个专用名词。 x(t)截断、周期延拓xT(t)周期信号xT(t)的傅里叶变换：对周期信号xT(t)采样，将离散序列xT(n),将积分转为集合：傅里叶变换公式按上式，用计算机编程很容易计算出指定频率点的值： f=？ /计算的频率点Fs=5120N=1024dt=1.0/Fspi=3.1415926 XR=0XI=0 For n=0 To N-1 XR=XR+x(n)*cos(2*pi*f*n*dt)*dt XI=XI+x(n)*sin(2*pi*f*n*dt)*dtNextA=sqr(XR*XR

7、+XI*XI) Q=atn(XI/XR)VBScript 样例展开，得连续傅立叶变换计算公式： 连续傅立叶变换编程计算实验： 采样信号频谱是一个连续频谱，不可能计算出所有频率点值，X( f )只能离散取值，频率取样间隔定义为：f=fs/N 频率取样点为0, f, 2f, 3f, .，有： 该公式就是离散傅立叶计算公式(DFT)2) 快速傅立叶变换 (FFT) FFT是实施DFT的一种有效算法，通过仔细选择和重新排列中间结果，速度上较DFT有明显的优点。展开各点的DFT计算公式：XR(1)=x(0).cos(2pi*0*1/N)+x(1).cos(2pi*1*1/N)+x(2).cos(2pi*

8、2*1/N).XR(2)=x(0).cos(2pi*0*2/N)+x(1).cos(2pi*1*2/N)+x(2).cos(2pi*2*2 /N).其中有大量cos、sin项的重复计算，FFT就是用技巧减少这些重复计算。 当采样点数为1024点,DFT要求一百万次以上计算量，而FFT则只要求10240次。 FFT的Matlab实现fs=1000 t=0:1/fs:0.6; f1=100; f2=300; x=sin(2*pi*f1*t)+sin(2*pi*f2*t); subplot(711) plot(x); title(f1(100Hz)f2(300Hz)的正弦信号，初相0) xlabel

9、(序列（n）) grid on number=512 y=fft(x,number); k=0:length(y)-1; f=fs*k/length(y); subplot(713) plot(f,abs(y); title(f1f2的正弦信号的FFT（512点）) xlabel(频率Hz) grid onx=x+randn(1,length(x); subplot(715) plot(x); title(原f1f2的正弦信号(含随机噪声) ) xlabel(序列（n）) grid ony=fft(x,number); k=0:length(y)-1; f=fs*k/length(y); su

10、bplot(717) plot(f,abs(y); title(原f1f2的正弦信号(含随机噪声)的FFT(512点) xlabel(频率Hz) grid on6 栅栏效应与窗函数 1)栅栏效应 采样信号的频谱，常采用FFT算法进行计算，设数据点数为： N = T/dt = T. fs 则计算得到的离散频率点为: Xs(fi) ， fi = k.fs/N , k = 0，1，2，.，N/2 相当于透过栅栏观赏风景，只能看到频谱的一部分，而其它频率点看不见，称此种现象为栅栏效应。 如果信号中的频率分量与频率取样点不重合，则只能按四舍五入的原则，取相邻的频率取样点谱线值代替。 栅栏效应误差实验：

11、2) 能量泄漏与栅栏效应的关系 频谱的离散造成了栅栏效应，谱峰越尖锐，产生误差的可能性就越大。 例如，余弦信号的频谱为线谱。信号频率与频谱离散取样点不等时，栅栏效应的误差为无穷大。 因信号截断产生了能量泄漏，即使信号频率与频谱离散取样点不等，也能得到一个该频率分量的近似值，因此能量泄漏误差不完全有害的。如没有信号时域截断产生的能量泄漏误差，频谱离散取样造成的栅栏效应误差不能接受。 能量泄漏分主瓣泄漏和旁瓣泄漏，主瓣泄漏可以减小因栅栏效应带来的谱峰幅值估计误差，有其好的一面，而旁瓣泄漏则是完全有害的。3) 常用的窗函数 (1) 矩形窗 (2) 三角窗 (3) 汉宁窗常用窗函数窗函数在减小栅栏效应

12、误差中的作用实验： 总结：信号截断能量泄漏FFT栅栏效应从克服栅栏效应误差的角度看，能量泄漏是有利的。通过加窗加大能量泄漏，减小栅栏效应误差：加矩形窗加BlackMan窗动手做：设计一个离散傅立叶计算程序，计算方波的频谱。观察栅栏效应带来的计算误差。7 常用的数字信号处理算法 数字信号处理是测试技术中最常用和最需要应用工程师掌握的部分，无论开发简单或复杂的测控系统或分析仪器，都会用到数字信号处理知识。 信号时域波形分析是最常用的信号分析手段，用示波器、万用表等普通仪器显示信号波形就可以特征参数。 1)时域波形参数计算 0At峰值/双峰值均值有效值方差周期波形分析的应用超门限报警 信号类型识别

13、信号基本参数识别 Pp-p 峰值P，双峰值Pp-ptAPPp-pP1=data(0)P2=data(0)For K = 0 To N If P1data(k) Then P2=data(k) End IfNextP=P1P2P=P1-P2 均值0AtU=0For K = 0 To N U=U+data(k)NextU=U/N 均方值E2=0For K = 0 To NE2=E2+data(k)*data(k)NextRMS=sqr(E2/N)方差V2=0For K = 0 To N V2=V2+(data(k)-U)*(data(k)-U)NextV=V2/N大方差 小方差 周期TtAT 穿越

14、判断法ATn=0AT=0.8*PFor K = 2 To N If data(k-1)AT And data(k-2)AT And data(k+2)AT Then ti(n)=K n=n+1 End IfNextT=(ti(2)-ti(1)*dt 两个序列的卷积：c(n)=a(n)*b(n)= 如a和b的脚标范围为：a1,a2与b1,b2，则c的脚标范围为a1+b1,a2+ b2 ，但matlab的索引是从1开始的，设a和b都是从n=1开始, 知道最后的范围即可.function c=convolution(a,b)M=length(a); %a序列的长度N=length(b); %b序列的

15、长度 c=zeros(1,N+M-1); %用来存放卷积的结果for n=2:(M+N) %卷积之后的索引本应该的范围 tmp_max=min(n-1,M); %根据求和表达式，由1=k=M 且 1=n-k x = 1:9; y=2:8; conv(x,y)ans = 2 7 16 30 50 77 112 147 182 197 200 190 166 127 72 convolution(x,y)ans = 2 7 16 30 50 77 112 147 182 197 200 190 166 127 72(3) 数字相关函数计算 变量之间的相依关系称为相关。信号之间的相似关系称为相关函数

16、。xyxyxy0At相关分析的工程应用 案例：机械加工表面粗糙度自相关分析 被测工件相关分析提取出回转误差等周期性的故障源。案例：地下输油管道漏损位置的探测 案例： AGV小车定位，声位笔定位实验：自相关分析如果a和b的脚标范围为：与，则c的脚标范围为，但是由于matlab的索引是从1开始的，这里假设a和b都是从n=1开始。知道最后的范围即可：function c =correlation(a,b)M = length(a);N = length(b);c=zeros(1,M+N-1); %用来存放卷积的结果for n=(-N+1):M-1 %相关运算本应该的范围 tmp_min = max(

17、1,1-n); %根据求和表达式，由1=k+n=M 且 1=k xcorr(x,y)=-0.0 -0.0 8.0 23.0 44.0 70.0 100.0 133.0 168.0 203.0 238.0 193.0 150.0 110.0 74.0 43.0 18.0 correlation(x,y) =8 23 44 70 100 133 168 203 238 193 150 110 74 43 18计算公式：dt=0.1;t=0:dt:100;x=cos(t);a,b=xcorr(x,unbiased);plot(b*dt,a) 对于互相关函数，把a,b=xcorr(x,unbiased

18、);改为 a,b=xcorr(x,y,unbiased);便可。 在matlab中实现自相关的无偏估计延迟序数(4) 功率谱估计与Matlab功率谱：随机信号的功率谱反映的是随机信号的频率成分及各成分的相对强弱。功率谱估计：基于有限的数据寻找信号、随机过程或系统的频率成分。分参数化和非参数化2大类方法。非参数化方法又称经典谱估计，根据截取的N个样本数据估计出其功率谱如周期图法、自相关法等，主要缺陷是描述功率谱波动的方差性能较差，频率分辨率低;参数化谱估计又称现代谱估计，如AR模型法、移动平均模型法(MA模型法)、自回归移动平均模型法(ARMA模型法)等 根据维纳辛钦定理，对于离散随机信号有:

19、Rx(m)为离散随机信号自相关函数，Sx(ej)为功率谱密度如获得Rx(m) ，Sx(ej)就不难求出, 故该法的实质是相关函数估值的傅立叶变换即为功率谱密度平稳随机信号的相关函数是确定性函数，其功率谱也是确定的 由平稳随机离散信号的有限个离散值x(0)，x(1)，x(N1)求出自相关函数然后在(-M，M)内对Rx(m)作Fourier变换，得到功率谱 用相关函数获得功率谱(BT法)功率谱可用FFT算法求出M趋于N时，上式求得的Rx(m)方差较大，可靠性差，谱估计质量下降为此Rx(m)上加窗函数w(m) 使方差大的自相关函数估值加权小，以减少对谱估计的影响功率谱是非负值，加窗后的功率谱也是非负

20、值因此窗函数的傅立叶变换必须是非负值，这样窗函数必须是偶对称的，上式中的三角形窗是常被采用 该法的缺点是M趋于N时，R(m)的方差很大，使谱估计质量下降；由R(m)得到的S(w)不一定为正值，从而可能失去功率谱的物理意义该法对由序列x(n)估计出的自相关函数R(m)进行傅立叶变换，得到x(n)的功率谱估计当延迟与数据长度相比很小时，有良好的估计精度Fs=500; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+0.1*randn(size(n); window=boxcar(length(xn);nfft=512;cxn=xcorr(xn, unb

21、iased);CXk=fft(cxn, nfft);Pxx=abs(CXk);index=0:round(nfft/21);k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1);figure(3);plot(k, plot_Pxx);基于Matlab实现的程序: Schuster首先提出的根据各态历经随机过程功率谱的定义进行的谱估计, 把随机序列x(n)的N个观测数据视为能量有限的序列，直接计算x(n)的DFT，得x(k)，然后取幅值的平方并除以N，作为序列x(n)真实功率谱的估计由平稳随机信号X(n)的有限个观察值x(0)，x(1)，x(N-1)求其

22、傅立叶变换然后进行谱估计 ，M=N1时，相关函数法与周期图法估计出的功率谱相同；MN1时，相关函数法的偏差大于周期图法,窗函数满足一定条件时是渐进无偏估计；其方差小于周期图法；分辨率比周期图法低，与窗函数的选择有关可见，窗函数的选择对谱估计的修正程度有影响周期图法(Periodogram) 周期图法因与序列的频谱有直接对应关系，且可采用FFT来计算而被广泛应用但需对无限长的平稳随机序列进行截断，相当于加矩形窗而成为有限长数据，这也意味着对自相关函数加三角窗，使功率谱与窗函数卷积，产生频谱泄漏，容易使弱信号的主瓣被强信号的旁瓣所淹没，造成频谱的模糊和失真，较低谱分辨率Nfft=512Nfft=1

23、024Fs=600;n=01/Fs1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+01*randn(size(n); window=boxcar(length(xn);nfft=512; Pxx，f=periodogram(xn，window，nfft，Fs);figure(1); plot(f，10*log10(Pxx);window=boxcar(length(xn); nfft=1024; Pxx，f=periodogram(xn，window，nfft，Fs);figure(2); plot(f，10+log10(Pxx); 周期图法基于Matlab实现的程

24、序: 周期图法不满足一致估计条件，需改进。改进目的主要是将周期图进行平滑，减小估计的协方差，得到一致谱估计。数据长度N太大时，周期图的功率谱估计的谱曲线起伏加剧；若N太小，谱的分辨率又不好。改进的主要方法有 Bartlett平均周期图法 将N点的有限长序列x(n)分段求周期图再平均将长度为N的数据分为L段，每段长度为M先对每段数据用周期图法进行谱估计，然后对L段求平均得到长度为N的数据的功率谱可得功率谱为 该估计的分辨率 式中由于M远小于N，故其分辨率比周期图法低可见，该法方差的改善是以牺牲分辨率为代价的 经典谱估计的改进Welch平滑平均周期图法 其改进了Bartlett法方差性能和分辨率基

25、本原理是先对随机序列分段，使每一段有部分重叠，然后对每一段数据用一个合适的窗函数进行平滑处理，最后对各段谱求平均得功率谱其中 ，w(n)是窗函数，其为矩形窗、重叠50%时，该估计的分辨率 Welch法允许各段数据交叠，数据段数L会增加，使方差得到更大的改善，但数据的交叠增加了各段数据的相关性，各段数据不再是完全独立的了,方差的减小不会达到理论程度另外采用合适的窗函数可以减少信号的频谱泄露，同时也可以增加谱峰的宽度，从而提高分辨率Fs=600;n=01/Fs1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n); nfft=512; window=boxcar(100); %矩形窗