用SP061A实现心电数据的FFT与压缩

在SP061A 单片机上实现对ECG信号的FFT、滤波和压缩。合理组织SP061A的硬件资源，并采取数据分段长度可选、避开高频分量的计算和简易的数据压缩算法，使存储开销、运算速度和精度满足实用要求。

在远程心电监护系统中，心电信号采集器是实现心电信号的现场采集、存储和传输的重要终端设备。对采集器的基本要求之一是：及时对采集到的心电信号进行滤波和压缩等预处理，以减少存储器占用量和数据远程传输到头端服务器的开销。为降低成本，这些任务一般采用单片机完成。然而，限于单片机的资源、运算能力和运行速度，许多压缩算法，如周期压缩法、小波变换压缩法和神经网络方法等无法使用，一些缺乏快速算法的频域变换法也很难达到实用的程度[3]。高性价比的心电信号采集器的研制一直是一个热点问题。
   
通过研究FFT（快速傅立叶变换）的算法结构和心电信号的特点发现，采用分段FFT，保留分析心电波形需要的谐波成分，巧妙地组织单片机的片内RAM资源，可使数据运算量和RAM开销大大减少，能实现数据滤波和压缩，且能达到实时采集与处理所需的运算速度。

SP061A 是凌阳科技公司研制的一款16位超低功耗单片机[1]，片内有2K字RAM、10位A/D转换器，CPU时钟高达49.152MHz，且价格低廉，还特别具有一套精简、高效的指令系统和类似于DSP的硬件内积运算功能。这些特点很适合心电信号的采集和处理。图1是作者研发的心电信号采集器中有关硬件的组成框图：多路ECG模拟信号送SP061A进行A/D转换，转换数据送NVRAM  DS1265W暂存；待采集完成后，由SP061A进行FFT和滤波、压缩；压缩结果送回DS1265W，再适时通过电话线或计算机网络送到监护中心处理、诊断。

本文仅讨论用SP061A实现FFT[2]、低通滤波与压缩。设对心电信号的采样率为500次／秒，数据精度为10位。

1  数据分段算法
   
设采集到的原始数据存于片外RAM中，将这些数据分为若干段，逐段读入片内进行FFT。各段的变换结果及时送回片外RAM中保存。
   
按照FFT的要求，段中包含的数据个数必须为2N，N为FFT变换的层数。考虑到SP061A片内RAM为2K字，此处取N=9或N=10，即段中数据为512或1024，以保证RAM够用。显然，段头和段尾的数据大小相等时，以该段作为一个周期而无限重复的波形将无跳跃点。经过“FFT变换到频域”→“丢弃高频成分”→“IFFT（快速傅立叶反变换，在头端PC上进行）”一系列操作而重建的时域波形，段与段之间的结合点将是连续的。但实际上，按上述分段几乎不能做到段头和段尾的数据大小相等。取两种段长的目的就是提供两种可能的选择——选择首尾数据之差较小的段作FFT。尽管如此，段首尾数据之差仍存在，经处理、复原后的波形在段的结合部位仍将有间断点。而采用加窗、延拓等办法在单片机上又难以实现。解决问题的策略为：分段时，各段间的数据首、尾各覆盖10个数据。头端PC在完成重建后，应将首、尾各5个数据丢弃。

2  时域数据的整序与加载
   
分段后，将该段加载到SP061A的RAM中，以实施FFT。原始数据以采集的时间先后顺序存放，加载时则应“整序”，即改变数据的先后顺序，以保证变换后的频域数据为正序。
   
设RS为指向片外RAM的、待加载的段内数据的偏移地址，RS=0…2N-1；Rd为指向片内RAM的、待写入数据的偏移地址，如图2。将RS按N位二进制逐位高低互换就得到Rd。例如，当N=9时，若RS为011001011B，则Rd为110100110B。为加快计算速度，将N=9时Rd的值制表存于FLASH ROM，供整序时查询。当N=10时，取RS的B0～B9位查表获得Rd，再将RS的B10位传送到Rd的B15位，最后将Rd循环左移1位。

FFT变换是复数运算。在将原始数据加载到片内RAM的同时，应把实数转换为复数，即令虚部为0。于是，一个原始数据加载到RAM中要占用2个字。复数的存储格式为：实部字存于低地址，虚部字存于相邻的高地址。现在考察RAM需要量。N=9时，段长为512个数据，加载到RAM中要占用 512×2=1024字；N=10时，段长为1024个数据，全部加载将占用 1024×2=2048字，超过片内RAM的可用容量。此时，将数据分为两部分，先将第一部分加载到RAM作FFT，得到中间结果，再将第二部分加载、变换，最后相加合成。

3  FFT变换及低通滤波
  
FFT将时域序列{x[i],i∈0…2N}变换为频域序列{F[i],i∈0…2N}。为了实现低通滤波，仅须保留{F[i]}中≤75Hz的频率分量。当N=9时，应保留{F[i]}中的前77个低频分量；当N=10时，则应保留{F[i]}中的前154个低频分量。这也同时减少了计算量，加快了计算速度；存放周转量所需的片内RAM也能得到保证。

为叙述简便，以N=3为例，研究FFT的计算结构，如图3所示。

kN=0…2N-1。将W[k]的实部和虚部都乘214，取整后制成表，存于FLASH ROM中，供程序查表获得其值；而W[k]与某数相乘，将32位运算结果右移14位作为积。这就使全部运算为整数运算，适应SP061A的硬件乘法功能。由图3知，第一层的计算仅涉及实部加减，虚部保持为0，可单独进行。从第二层开始有复数乘，但是，当只需计算{F[i]}中的低频分量时，许多中间结果可不计算。例如，如果需计算出F[0]和F[1]（即保留原始信号的直流分量和1次谐波），则仅需计算x[0]3、x[4]3和x[1]3、x[5]3。计算层数N越多，减少的运算也越多。
   
复数乘可利用SP061A的内积功能实现。例如，要计算x[i]×W[j]，设x[i]×W[j]=(a+jb)×(c+jd)=ac+(-bd)+j(bc+ad)。显然，结果的实部和虚部均为内积形式，只是设置操作数时须注意符号和排列顺序。
   
上述方法使计算量显著减少。以512点FFT为例，计算出全部频率分量需要512×log2512=4608次运算，其中含有2048次复数乘。若计算77个低频分量，则只有3611次运算，其中含有1767次复数乘。
   
当N=10时，计算点数达1024，片内RAM不够用。此时，应按1024点的整序次序取数，先对x[0]1～x[511]1进行FFT，算出F1[0]～F1[153]，暂存于片内RAM中的一个缓冲区；再对x[512]1～x[1023]1进行FFT，算出F2[0]～F2[153]；则最终结果为:F[i]=F1[i]+F2[i]，i=0…153。
   
为避免计算中产生数据溢出，从第三层开始，对x[i]4～x[i]9都算术右移1位。操作的累积结果使F[i]缩小了64倍，故在重建时应扩大64倍。如此操作实际上降低了运算精度，但实验表明，重建的波形完全满足医学观察要求。

4 数据压缩
   
采取如下简易格式实现数据压缩：
   
对于F[0]，因虚部为0，仅用一个字存放实部，重建时默认虚部为0；

对于F[i],i>0，若实部在-64～63范围内且虚部在-128～127范围内，则用2个字节存放，格式如下：

两种格式由第1字节的最高位区分。

5 实验结果与分析
   
用自行研发的心电信号采集器进行实验，对采集到的4个样本进行处理，实验结果如表1。表1中，PRD为均方根误差，CC为相关系数，计算公式为：

数据系列及其平均值，平均值。

处理时间为SP061A完成FFT与压缩花费的时间, CPU时钟设置为49.152MHz。
   
实验表明，本方法用价格低廉的单片机实现了复杂的FFT与数据压缩，计算耗时少，所得结果满足实用要求。由图4可见，重建后的波形在段间结合点无畸变。噪声较弱时PRD和CC参数较为理想；而当噪声很强时，如图4（a）、4（b），因滤除了高频噪声而使得重建波形与原始波形差距较大，PRD和CC参数已不能说明问题。压缩算法简便，CR约为4。顺便说明，本方法未实现50Hz干扰滤波、肌电干扰滤波和基线漂移，这些处理可在头端PC上进行。

参考文献

1 北阳电子内部技术资料.SPCE061A原理与应用教材书. http://www.unsp.com.cn. 2003.5
     
2 胡 阳， 杨叔子. 一种基于Fourier复变换的压缩比自调整ECG数据压缩算法[J]. 上海：中国医疗器械杂志,1997；21(2):75~78
      
3 张浙亮，吕维雪. 心电信号数据压缩技术的发展[J]. 国外医学:生物医学工程分册,1997；20(2)：73~80

4 王培康，费小英.利用DCT分量差值压缩ECG数据的方法[J]. 中国生物医学工程学报,2002；21(5):456~460