一种振动信号处理过程的数据量估计方法
摘 要：采用硬件板卡在对振动信号进行实时处理时，需要对信号处理过程的运算量进行估计，进而作为选择处理器和相关硬件配置的依据。本文从Cache使能、傅里叶变换、位翻转和读写片外SRAM等四个方面，精确估计了一种振动信号进行FFT运算时的数据量和时间。该方法的估计结果与实测数据一致，可作为相关算法运算量估计的参考。
关键字：振动信号；数据量估计；快速傅里叶变换
引言
安装在飞机或台架上的航空发动机，是一个无限多自由度的振动系统，发动机的整机振动就是这一系统在各种激振里作用下的响应[1]。整机振动是影响发动机寿命和飞行安全的决定性因素[2]，对发动机整机振动进行监控，可以为发动机正常且安全使用提供有力保障。
通常对振动信号处理的方法有：时域分析法、频域分析法及时频分析法。时域分析是采用相关分析法和特征值分析法，根据信号随时间变化的波形来分析信号，从而得出信号的特征量和组成。频域分析法是采用频谱分析法和倒谱分析法，将信号从时域转换到频域，分析频谱的组成，得到信号的本质特性，频域分析法虽能反映信号的本质，但不能体现出信号频率局部瞬间的变化，具有局限性。时频分析法通过短时傅里叶变换、小波变换等方法，是针对非线性非平稳信号提出的处理方法。
傅里叶变换是时域转换到频域的基本工具，是目前频谱分析的主要方法，而快速傅里叶变换（FFT）则是巧妙利用W因子的周期性和对称性推导得到的一个高效的离散傅里叶变换的快速算法[3]，在振动信号处理中被普遍使用。本文对振动信号从时域转换到频域使用的FFT运算，给出了该运算所需的数据量和时间估计方法，对系统使用的运算量进行了精确估计，为选择处理器或者硬件设计提供了依据。
2 振动信号的快速傅里叶变换原理
有限长离散信号x（n）（n=1，2，…，N-1）的离散傅里叶变换定义为：
可以看出，离散傅里叶变换需要计算大约N2次乘法和N2次加法。利用的WN对称性和周期性，将点离散傅里叶变换分解为两个N/2点离散傅里叶变换，总的计算量只是原来的一半，即
将x（n）分解为偶数和奇数的两个序列之和，即x（n）-x1（n）+x2（n）。x1（n）和x2（n）的长度都是N/2，x1（n）是偶数序列，x2（n）是奇数序列，则
所以：
由于
其中X1（k）和X2（k）分别为x1（n）和x2（n）的N/2点的离散傅里叶变换。由于X1（k）和X2（k）均以N/2为周期，且
上式的运算可以用图1表示，根据其形状称之为蝶形运算。依次类推，经过m-1次分解，最后N将点离散傅里叶变换分解为N/2个两点的离散傅里叶变换。
图1 蝶形运算
快速傅里叶变换的原理就是通过许多小的更加容易进行的变换去实现大规模的变换，降低了运算要求，提高了运算速度。
3 振动信号处理的软硬件实现
3.1 硬件结构
本系统的硬件平台主要由PowerPC处理器、DSP处理器、FPGA，以及相关的外围电路组成，系统原理参见图2。
图2 系统原理图
本系统完成的主要功能是，FPGA控制AD转换器对振动信号进行采集，数据在其内RAM中缓存后被发送至DSP，DSP对振动信号进行傅里叶变换处理，获取频域信息，并对信号特征（加速度、速度及位移）进行提取，之后将振动信号的频谱数据和信号特征送给FGPA，FPGA通过双口RAM再将相关数据实时传送给上位机，供用户使用。
3.2 信号处理流程
本系统中的DSP软件完成的功能如下：
a.从FPGA的FIFO中获取振动信号数据，存储在片外RAM中；
b.产生傅里叶变换所需要的系数，存储在内存中；
c.从片外RAM中提取16384个数据，进行快速傅里叶变换；
d.将快速傅里叶变换的结果，通过DPRAM传送出去。
振动信号处理流程参见图3。
图3 振动数据处理流程图
4 信号处理中的运算量分析
4.1 理论运算量分析
根据图3的处理流程可知，第2、3和4个步骤构成单路数据快速傅里叶变换。根据流程分析，可从Cache使能的影响、傅里叶变换所需执行时间、读写片外SRAM需要的时间及结果位翻转的时间等几个方面，对所需要理论时间进行估计。
4.1.1 Cache使能的影响
DSP6713内部有两级的程序和数据存储器，L1P为一级高速程序缓冲存储器、L1D为一级高速数据缓冲存储器。L2为二级存储器，二级的程序和数据存储器可共用。
对于数据的访问策略：
当CPU访问的数据未命中时，L1D会将请求的数据替换一个新的行，若L1D未命中但命中了L2 SRAM或L2，高速缓存器会给CPU带来4个周期的延迟（假设没有其他的L2访问延迟L1D的数据访问）；若L1D未命中并且没有命中L2时，会从外部存储器中读取数据，一旦这些数据读回，它将被存储在L2中并传输到L1D，外部访问未命中带来的不利因素根据外存宽度的不同和其他系统负载的不同而不同。
对于程序的访问策略：
命中L1P的程序取指，在一个周期内完成，不会给CPU带来延迟。在没有其他存储器操作对L2访问产生延迟的情况下，不管相关代码是否并行，L1P未命中而L2命中会给CPU带来5个周期延迟。在L1P和L2都未命中的情况下，产生的延迟会持续到L2从外部存储器读回程序并传送到L1P。产生的延迟依赖于外部程序存储器的不同和系统负载的不同而不同。
Cache使能的状态下，当CPU发出访存请求后，存储器地址先被送到Cache控制器以确定所需数据是否已在Cache中，若命中则直接对Cache进行访问。
本课题Cache设置为：
*CCFG=0x07；
*MAR0=0x01；
CCFG被配置为0x07，即使能64K大小的L2Cache，MAR0被配置成0x1，即相应的存储器范围可Cache。
4.1.2 傅里叶变换所需执行时间
快速傅里叶变换所需的时间为复数乘法和复数加法所耗费时间之和，根据第二章的分析可知，对于N个点的傅里叶变换所需的复数乘法数为（N/2）log2N，复数加法数为Nlog2N。
对于N=16384，傅里叶变换完全在内存中进行，根据计算可知，需要复数乘法运算114688次，复数加法运算229376次。
在《DSP算法与体系结构实现技术》书中提到，为实现流水线的正常工作，乘法累加运算往往需要在一个时钟周期内完成，一次复数乘法需要4次乘累加，故针对N为16384个点的傅里叶变换，所需要的时间T的计算公式如下：
T=乘累加次数*执行时间=Mc*4*1fs （7）
式中：Mc为复数乘法数，fs为采样频率。
由公式（7）可得到，若内存足够大，16384个点的傅里叶变换需要耗时2.29376ms。
4.1.3 读写片外SRAM需要的时间
对于DSP6713，片外SRAM的访问，是通过外部的存储器接口（EMIF）完成，因此EMIF的速率决定了，片外SRAM的访问速率。因此，首先要配置锁相环PLL，来产生供EMIF使用的时钟；其次，配置EMIF，设定读写SRAM的时间；最后根据程序的调用，计算读写SRAM的时间。
a.设置PLL控制器
本课题PLL设置为：
st_Pll.pllm = 0xa；
st_Pll.plldiv0 = 0x8001；
st_Pll.plldiv3 = 0x8003；
st_Pll.pllcsr = 0x0；
系统输入频率为40MHz，通过PLL设置后，产生过程如下：首先，经过D0分频器，分频系数为2，分成20MHz；其次，经过PLLM倍频器，倍频系数为0xa，倍频为200MHz；最后，经过D3分频器，分频系数为4，分成50MHz，供EMIF使用。
b.配置EMIF
在本课题中，EMIF的配置如下：
*（int *）EMIF_GCTL = 0x00003078；
*（int *）EMIF_CE0 = 0x0081C521；
*（int *）EMIF_CE1 = 0x1151c511；
*（int *）EMIF_CE2 = 0x10914611；
*（int *）EMIF_CE3 = 0x10914521；
*（int *）EMIF_SDRAMCTL = 0x5748F000；
*（int *）EMIF_SDRAMTIMING = 0x005DC5DC；
*（int *）EMIF_SDRAMEXT = 0x0017DF3F；
根据以上的配置可得到：写建立宽度为1个时钟周期；片外SRAM为32位宽的异步接口；选择32位宽的异步接口；读选通宽度为5个时钟周期；最小转换时间为3个时钟周期；读建立宽度为1个时钟周期；写保持宽度为0个时钟周期；写选通宽度为2个时钟周期；写建立宽度为0个时钟周期。
c.计算读写SRAM的时间
如图4所示，调用DSPF_sp_cfftr2_dit函数执行快速傅里叶变换过程需执行循环14次，将16384*2个点原始点（实部为原始点，虚部为0）和对应的系数16384*2（复数），送入DSPF_sp_cfftr2_dit函数进行运算，并将计算的结果16384*2（实部+虚部）存入片外SRAM中。读选通需要5个EMIF时钟，写选通需要2个EMIF时钟，因此需要时间如下：
原始点读入时间：
16384*2*14*5/50MHz = 45.8752ms
旋转因子读入时间：
16384*2*14*5/50MHz= 45.8752ms
FFT结果读出时间：
16384*2*14*2/50MHz = 18.3501ms
图4 FFT库函数的C实现
4.1.4 结果位翻转需要的时间
图5为位翻转函数，在整个快速傅里叶变换过程中，需要包含两次位翻转：
a.在产生了系数之后，按照蝶形运算，需要将系数进行位翻转；
b.在调用库函数，完成FFT变换之后，需要将逆序的结果变为正序。
图5 位翻转函数
本函数需要循环16384次，读选通需要5个EMIF时钟，写选通需要2个EMIF时钟，系数位翻转及FFT结果位翻转均需要从外部将16384*2（包括实部和虚部）个系数或结果，从外部SRAM中搬移到Cache及内存中处理，处理完成后再写入外部SRAM中，因此需要的时间如下：
系数位翻转时间：
16384*2*5/50MHz + 16384*2*2/50MHz = 4.5875ms
结果位翻转时间：
16384*2*5/50MHz + 16384*2*2/50MHz = 4.5875ms
4.15 计算过程
FFT运算需要时间 = FFT理论执行时间 +原始点读入时间 +旋转因子读入时间 + FFT结果读出时间 + 系数位翻转时间 + 结果位翻转时间
= 2.29376 + 45.8752 + 45.8752 + 18.3501 + 4.5875 + 4.5875
= 121.5693ms
4.2 实际运算量分析
按照图3的方法设计软件代码，在图2的硬件环境中运行，使用示波器测量，可得到：单路数据进行快速傅里叶变换，包括调用库函数，系数产生及翻转，总共需要128ms，基本与估计值一致。
5 结论
本文所描述的振动信号处理过程的数据量估计方法，基本可估计出整个过程中，数据量和资源的使用情况，可作为选择处理器或者硬件设计的基础。
除此之外，若想将运算速度提高，可从以下几个方面着手改进：
1.从软件角度，将EMIF的速度从50MHz提升到100MHz，可以提高约50%的速度；
2.可以将变换系数在初始化阶段生成，并作为固定数据表格放入内存中，约占128KB，使用时，查表获取即可。
3.可以采用速度更高或内存更大的DSP作为CPU，进行大数据量的FFT变换的处理器。
参考文献：
[1] 《航空发动机设计手册》编委会. 转子动力学及整机振动[K]. 北京：航空工业出版社，2000.
[2] 李宝凤，王德友. 某型航空发动机300h持久试车中的整机振动分析[J]. 航空发动机，2004. 19-21.
[3] Cooley J W， Tukey J W. An Algorithm for the Machine Calculation of Complex Fourier Series[J]. Mathematics of Computation， 1956，19（90）：297-301.
[4] 许邦建，孙永节，唐涛. DSP算法与体系结构实现技术[K]. 北京：国防工业出版社，2010.
作者简介：张蕾，女，1984年12月，硕士，工程师，健康管理方向
