频谱分析中数据处理的几个问题

2008 AVL先进模拟技术中国用户大会论文 1 图 图图 图1 离散 离散离散 离散正弦时间序列 正弦时间序列正弦时间序列 正弦时间序列信号 信号信号 信号(总长 总长总长 总长100点 点点 点) Time History -10 -5 0 5 10 Si gn a l 0 0.2 0.4 0.6 0.8 1 Time [s] 频谱中数据处理的几个问题 蓝军 （重庆长安汽研院） 摘 要：讨论了 NVH数据处理中常考虑的几个重要问题。 关键词：NVH，FFT，结果数据步长，采样频率，频率分辨率，信号分析频率 主要软件：AVL EXCITE，Impress Chart，MS Excel 1. 前言 通常 NVH动力学计算有大量的数据需要处理，用来分析激励和响应的信号特征，结合 模态分析，寻找减振降噪的对策。在计算机出现后，相对于模拟信号处理，现代数字计算方 法得到了长足发展，其中经常用到离散信号的频谱分析、相关分析、滤波等。计算与试验数 据进行对比，也需要关注数据处理的方式和要求。 2. 离散信号 为校核数字算法是否适用于某种信号，在分析一个离散数据信号时，经常需要产生调试 用的一个正弦时间序列信号（因为一个信号总可以展开成多个正弦信号叠加）： x(n*ts)=A*sin(2*pi *fo*ts*n+ϕ) (1) 其中:n=1,2,…N，（一般希望数据采样点数总长为 N=2k，k为整数，比如：N=1024点）， A为信号幅值， fo为正弦信号频率， ts为采样周期（秒），就是在计算中的 结果存盘步长（秒为单位）， 因此，采样频率 fs=1/ts，单位 Hz， 频率分辨率 ∆fanalysis=fs/N， 常用 EXCEL来产生正弦信号的离散时 间序列（pi用 PI()示），图 1为用 fs=100Hz生成的一个频率是 fo=3Hz幅值为 10的正弦时 间序列。为视觉着想，常用连线表示离散时间序列。 3. EXCITE的结果存贮步长和信号最高分析频率 根据采样定理，采样频率 fs ≥2×信号最高分析频率 fanalysis，一般可取为 2.5×fanalysis， 比如要分析 3000Hz以内的噪声信号，最低需取采样频率 fs＝7500Hz。 下面推导采样时间步长 ts与发动机曲柄转角步长 ta之间的关系： � 对于发动机转速为 n(rpm)，旋转一周的时间＝60/n (s/cycle)， � 使用 fs(Hz)，其采样周期 ts=1/fs (s/一次采样)， � 则发动机旋转一周的采样次数＝(60/n)/(1/fs)=60*fs/n (cycle/一次采样)， 因此，数据存盘的曲柄转角步长 ta＝360/(60*fs/n)=6*n/ fs (CAdeg/cycle) (2) 在 EXCITE的 simulation control中仿真步长一般取得都较小（如 0.0625deg），用来满足 迭代收敛要求。需注意：如果 simulation control->data storage和 result control中结果存盘步 长 ta＝1deg，则可处理的信号最大分析频率与发动机转速的关系如下表： 2008 AVL先进模拟技术中国用户大会论文 2 图 图图 图2 EXCITE PU设定的数据长度对频率分辨率 设定的数据长度对频率分辨率设定的数据长度对频率分辨率 设定的数据长度对频率分辨率的影响 的影响的影响 的影响 -0.003 0 0.003 0.006 0.009 0.012 0.015 0.018 D is pl a ce m e n t[m m ] 0.02 0.024 0.028 0.032 0.036 0.04 Time[s] FFT 0 0.002 0.004 0.006 0.008 0.01 0.012 M a gn itu de [m m ] 0 200 400 600 800 1000 Frequency[Hz] -0.003 0 0.003 0.006 0.009 0.012 0.015 0.018 D is pl a ce m e n t[m m ] 0.02 0.04 0.06 0.08 0.1 0.12 Time[s] FFT 0 0.002 0.004 0.006 0.008 0.01 0.012 M a gn itu de [m m ] 0 200 400 600 800 1000 Frequency [Hz] 6000rpm一个循环结果 一个循环结果一个循环结果 一个循环结果 6000rpm五个循环结果 五个循环结果五个循环结果 五个循环结果 发动机转速（rpm） 1000 2000 3000 4000 5000 6000 信号最高分析频率 fanalysis（Hz） 2400 4800 7200 9600 12000 14400 进行 EXCITE的振动噪声计算，2000rpm以上转速的结果存盘步长 ta=1deg，则可以满 足最高分析频率的要求。 4. FFT的数据长度和频率分辨率 FFT算法是 Cooley 和 Tukey于 1965年，巧妙编制的一种计算离散傅立叶变换 DFT的 快速算法，Impress Chart的运算中(Operations)提供了快速傅立叶变换 FFT，可获取幅值谱和 相位谱（没有提供功率谱密度，可自己处理产生）；也提供了快速傅立叶逆变换运算 Inverse FFT。FFT得到的频谱是离散的（图 3图 4），为视觉着想，图中可用连线表示。 如果希望细化了解信号的频率成分，可缩小频率分辨率 ∆fanalysis=fs/N，其中增加数据采 样点数是有效的办法。图 2 为一 6000rpm 的结果，使用一个循环结果时能识别的最小频率 分辨率是 50Hz，使用 5个循环时频率分别率可达到 10Hz。如果处理试验数据，并关注 0～ 200Hz范围内的振动问题时，可选用适当的采用频率和尽可能厂的数据长度。 1个循环结果 个循环结果个循环结果 个循环结果：存盘步长 1deg CA 5个循环结果 个循环结果个循环结果 个循环结果：存盘步长 1deg CA 采样点数总长 N=720 采样点数总长 N=720×5 采样频率 fs＝6×6000/1 (Hz) 采样频率 fs＝6×6000/1 (Hz) 频率分辨率∆f=50Hz 频率分辨率∆f=10Hz 2008 AVL先进模拟技术中国用户大会论文 3 Time History -10 -5 0 5 10 Si gn a l 0 0.2 0.4 0.6 0.8 1 Time [s] 整整整 非整整整 FFT 0 2 4 6 8 10 M a gn itu de 0 3 6 9 12 15 Frequency [Hz] 整整整 非整整整 图 图图 图3 数据不完成周期的影响 数据不完成周期的影响数据不完成周期的影响 数据不完成周期的影响 Time History (补 补补 补补 补补 补补 补补 补补 补补 补) -10 -5 0 5 10 Si gn a l 0 0.2 0.4 0.6 0.8 1 Time[s] 不不整整整 补补补补 FFT 0 2 4 6 8 10 M a gn itu de 0 3 6 9 12 15 Frequency [Hz] 不不整整整 补补补补 图 图图 图4 数据补零的影响 数据补零的影响数据补零的影响 数据补零的影响 图 2振动位移曲线的 FFT分析还反映出：在发动机 6000rpm下，基频 f 基 = 6000rpm/60s = 100Hz，四冲程发动机存在 0.5谐次主频（对应于 50Hz），因此使用 1个循环结果正好可 以反应相应谐次的情况。但是如果希望详细了解发动机 6000rpm 下的凸轮轴承力的情况， 可在 TYCON中输出 5个循环结果，频谱的频率分辨率则可为 10Hz，则 10Hz的频率分辨率 能将其频率响应更清晰地表现出来。对于试验结果的处理，由于有干扰信号的影响，为提高 信噪比，也需要采用包含众多循环的数据来进行统计。 5. FFT的频谱旁瓣泄漏 FFT 算法要求信号数据截断是完整周期的（即要求数据长度正好为一个或几个信号周 期），这样离散时间序列才可精确地在整个时域内周期性地延拓，为此我们输出数据都采用 整个循环结果。图 3红色数据是图 1中的完整周期性的正弦时间序列，FFT频谱幅值为 10， 频率为 3Hz，蓝色数据为非完整周期的数据，经 FFT 分析泄漏到其它频率成分中，还偏移 了主频的主瓣（往频率减小或增大的移动都有可能，与数据截取状态有关）。图 4中表明， 将信号补零处理到完整周期长度后，虽然有所改善，但由于信号的缺损，效果也不理想。 上面图 2中，EXCITE计算信号无论是一个循环还是 5个循环，都是完整周期的，因此 FFT 计算的泄漏较小，因此进行 FFT 处理时常可不加窗。应避免使用非整循环的数据来进 行 FFT（比如让 EXCITE输出 5.5个循环的数据）。试验数据处理时，很难正好采样到整循 环，因此减小泄漏就显得更加重要，一般需加窗。 6. FFT的窗函数 FFT 分析中常常要用到窗函数，用来减小频谱泄漏。要清楚的是，“截断数据”本身就 是对原信号加矩形窗(也称为“Uniform窗”或“不加窗”)。大家知道，窗函数的频谱由一个 主瓣和几个旁瓣组成，旁瓣在主瓣的两侧以一定的间隔衰减至零。如果时间序列的长度不是 周期的整数倍，窗函数的连续频谱将偏离主瓣的中心，频率偏移量对应着信号频率和 FFT 频 率分辨率的差异，这个偏移导致了频谱中出现旁瓣，所以，窗函数的旁瓣特性直接影响着各 2008 AVL先进模拟技术中国用户大会论文 4 图 图图 图5 泄漏和加窗 泄漏和加窗泄漏和加窗 泄漏和加窗 图 图图 图6 窗函数对窄带信号 窗函数对窄带信号窗函数对窄带信号 窗函数对窄带信号的影响 的影响的影响 的影响 频谱分量向相邻频谱的泄漏宽度。 表 1和表 2是常用窗函数特性及其选用要求。要选择合适的窗函数，须先估计信号的频 谱成份，如果信号含有频率很接近的两个或两个以上成分，应选用主瓣很窄的矩形窗；如果 强干扰频率分量紧邻被测频率时，应选择旁瓣峰值较小的窗函数。 表 1 常用的窗函数特性 窗类型 -3dB主瓣宽度 -6dB主瓣宽度 旁瓣峰值 旁瓣衰减速度 矩形窗 Uniform 0.89 bin 1.21 bin -13 dB 20dB/十倍频程 6dB/八倍频程 汉宁窗 Hanning 1.44 bin 2.00 bin -32 dB 60dB/十倍频程 18dB/八倍频程 表 2 信号频谱与窗函数的选择 窗类型 适用于信号类型 窗特性 矩形窗 Uniform 平稳宽带随机 分段逼近正弦波 优点：主瓣窄。缺点：旁瓣衰减速度低、频率 分辨率较低 汉宁窗 Hanning 平稳窄带随机 正弦波或正弦波组合 优点：瓣幅峰值较大、频率分辨率高、可减小 泄漏、旁瓣衰减速度快。缺点：主瓣稍宽 图 5中人工产生一个信号 x(n*ts)=10*sin(2*PI()*140*ts*n)+6*sin(1.3+ 2*PI() *180*ts*n)， 是两个正弦函数的叠加，其中一个幅值 10频率 140Hz，另一个幅值 6频率 180Hz，ts=0.0001s， 数据总长 1024 点。该采样信号 0.1024s，明显不是完整的周期（首尾相连的情况就可见）， 因此有较大泄漏。图中反映了不加窗和加 Hanning窗的幅值谱对比，可见不加窗的泄漏较为 严重，140Hz幅值仅 8.5，180Hz的幅值仅 4.5；加 Hanning窗后，泄漏减小，140Hz幅值仅 9.5，180Hz的幅值仅 5.4，更接近真值。 图 6为 EXCITE PU计算得到的曲轴飞轮端横向位移（5个循环），虽然呈周期性，但它 具有平稳窄带随机特征，且频率成分相距很近，为减小频率成分间的干扰，可使用 Hanning 窗泄漏小，各频率的幅值更加精确。同时也可获得良好的坎贝尔图，如图 7和图 8。 泄漏大 泄漏小 2008 AVL先进模拟技术中国用户大会论文 5 图 图图 图7 坎贝尔图 坎贝尔图坎贝尔图 坎贝尔图 图 图图 图8 FFT时的起止时刻 时的起止时刻时的起止时刻 时的起止时刻 图 9上图为 EXCITE TD计算得到 的气门落座力曲线（2个循环结果），可 见它具有宽带随机特征，一个明显的特 征就是该信号含有众多的幅值相对很 大的脉冲。气门落座冲击力越大，激励 出的宽频噪声频率成分也越丰富，使缸 头辐射噪声的强度、频率范围和声音品 质都恶化。中图和下图是幅值谱，为清 楚说明问题，采用了两种作图方式，正 确的理解应该是下图的离散幅值谱。图 中可以看到加 Hanning 窗后，由于 Hanning 窗的主瓣过宽，相邻频率的成 分互相干扰，致使弱信号成分被淹没 了。关于是否加窗，一定要先了解时域 信号特征，再决定加窗类型。 7. FFT的处理增益 EXCITE计算数据中不存在干扰信号，但如果要将 EXCITE频谱与试验频谱进行比对， 那么试验频谱的处理应该注意提高其增益，以便得到满意的比对效果。 图 10和图 11是用 EXCEL产生的调试数据，它就是图 1中的幅值为 10频率为 3Hz的 正弦信号叠加了随机干扰信号：x(n*ts)=10*sin(2*pi()*3*0.01*n)+3*rand()，其中，ts=0.01s。 图 10 中 n=1,2,…100，可见，随机干扰信号在 3Hz 以外的频率下产生干扰分量，经 5 次平 均后，有明显改善，图 11 中增加每段数据长度 n=1,2,…500，使 3Hz 以外的干扰降低到了 对数 0.1以下，显著地提高了处理增益，也提高了频率分辨率。设备处理时可使用 64次或 128次平均，要求信号是在稳定工况下平稳随机的，平均次数多了会增加试验的总时间，注 FFT 0 50 100 150 200 250 300 350 M a gn itu de 0 200 400 600 800 1000 Frequency [Hz] Hanning None -500 0 500 1000 1500 2000 2500 Ca m Sh a ft Be a rin g Fo rc e [N ] 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.055 0.06 time[s] z-Force[N] FFT 0 50 100 150 200 250 300 350 M a gn itu de 0 200 400 600 800 1000 Frequency [Hz] Hanning None 图 图图 图9 窗函数对宽带信号的影响 窗函数对宽带信号的影响窗函数对宽带信号的影响 窗函数对宽带信号的影响 脉冲 脉冲脉冲 脉冲 2008 AVL先进模拟技术中国用户大会论文 6 Sine with Random signal -10 -5 0 5 10 15 Si gn a l 0 0.2 0.4 0.6 0.8 1 Time[s] s1 s2 s3 s4 s5 FFT 0.01 0.1 1 10 M a gn itu de [对 对 ] 0 5 10 15 20 25 30 35 40 45 50 Frequency [Hz] s1 s2 s3 s4 s5 平平 Sine with Random signal -10 -5 0 5 10 15 Si gn a l 0 1 2 3 4 5 Time[s] s1 s2 s3 s4 s5 FFT 0.01 0.1 1 10 M a gn itu de [对 对 ] 0 5 10 15 20 25 30 35 40 45 50 Frequency [Hz] s1 s2 s3 s4 s5 平平 图 图图 图10 平均处理 平均处理平均处理 平均处理 图 图图 图11 平均处理 平均处理平均处理 平均处理and增加数据长度 增加数据长度增加数据长度 增加数据长度 图 图图 图12 倍频程和 倍频程和倍频程和 倍频程和1/3倍频程 倍频程倍频程 倍频程 1/3 Octave Bands 意协调。 8. 倍频程和 1/3 倍频程 为宏观地了解一定频率范围内的振动噪声状态，可使用倍频程或 1/3倍频程（图 12）。 Impress Chart中对原信号进行 FFT后（不用加窗函数），再转成相应倍频程谱图即可。对于 信号可采用分贝来表述其幅值大小： )0/log(20)( YYdBY = Y0为参考值，AVL常取：对速度 v0 = 5.e-8 m/s，对加速度 a0 = 1.e-6 m/s2 图 13是调试正弦函数 F(n*ts)=10*sin(2*PI()*n*f0*ts)，幅值为 10，固有频率 f0=600Hz， 采样周期 ts=0.0001s的时域、频域图，其中 600Hz分别对应第 17号 1/3倍频程带，或第 6 倍频程带。使用 Impress Chart 进行 FFT 变换时不需加窗函数，再使用 Frequency Bands， Method选用 Integral Levels，转换成倍频程或 1/3倍频程。 当取参考值 Y0=1，幅值 Y=10的正弦函数幅值分贝级=20*log(10/1)=20(dB)，图 13中可 2008 AVL先进模拟技术中国用户大会论文 7 -10 -5 0 5 10 A m pl itu de 0 0.02 0.04 0.06 0.08 0.1 Time[s] -50 -40 -30 -20 -10 0 10 20 M a gn itu de dB (20 * lo g( y/ 1) 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Frequency [Hz] FFT without window 图 图图 图13 调试信号的处理对比 调试信号的处理对比调试信号的处理对比 调试信号的处理对比 1/3 Octave 0 5 10 15 20 M a gn itu de dB (20 * lo g(y /1 ) 5 10 15 20 25 1/3 Octave Bands Octave 0 5 10 15 20 M ag n itu de dB (20 * lo g(y /1 ) 1 2 3 4 5 6 7 8 9 Octave Bands 见在 1/3倍频程中虽然有微量的泄漏，但处理误差并不大。 9. 结语结语结语结语 1、信号调试和处理是实践性非常强的技术，需要根据信号本身的特征来确定处理的参 数定义，才能获得准确而清晰的分析。 2、信号处理包括激励分析、响应分析、传函分析等，做为减振降噪信号分析手段。 参考文献 [1] (美)Richard G. Lyons著,朱光明译, "数字信号处理（原第 2版）", 机械工业出版社,2006年