1 引言
科里奥利质量流量计（以下简称科氏流量计）可直接高精度地测量流体的质量流量且可同时获取流体密度值，是当前发展最为迅速的流量计之一。
科氏流量计通过测量2个磁电传感器输出信号之间的相位差（时间差）来得到流体的质量流量。为了能够从含有噪声的传感器输出信号中准确地提取出流量信息，国内外研究者采用数字信号处理技术来计算相位差/时间差。但均是以时不变信号为研究对象的，即认为在一批处理数据（例如2048点）中，信号的频率、相位和幅值是不变的。实际上，由于受到管内流体的流速、密度和流体脉动等因素变化的影响，科氏流量计信号会随着时间发生变化。为此，我们率先建立了时变信号模型，并提出了基于陷波滤波器和滑动Goertzel算法的科氏流量计信号处理算法，克服了基于DFT频谱分析整周期采样的影响，在每个采样点都可输出计算结果，并且能跟踪变化的流量信号[1-2]。但是，算法的收敛过程较长，跟踪精度不高，且跟踪不到信号的微小变化，同时用于实际系统时计算量大且会发生数值溢出。在此基础上，文献[3-4]提出了计及负频率影响的DTFT递推算法，即在计算正弦信号的傅里叶系数时不忽略负频率成分的影响，从而提高计算精度，极大地缩短了收敛过程。但该算法只适合一段时间内信号恒定的情况，对于流量变化如批料流或有微小波动的时变信号时就失去了作用。同时，文献[2-3]没有给出算法的具体实施方案、DSP芯片的型号、信号处理硬件系统和采样频率等信息。
为此，针对时变信号，提出将多抽一滤波、自适应格型陷波滤波和负频率修正的滑动DTFT递推算法（以下简称SDTFT）有机地组合起来形成一套完整的科氏流量计信号处理算法，并在DSP系统上实现，实时处理流量计信号，既能用于非时变信号又能跟踪频率、相位变化的信号，更加接近实际应用情况。且算法计算量较小，不会发生数值溢出，便于实时实现。
2 时变信号
根据科氏流量计的工作原理，在理想情况下，传感器输出的信号为标准的正弦信号。但在实际应用中，由于受到管内流体特性如流速、密度、流体脉动及流场等因素变化的影响，信号的频率和相位并不是恒定不变的。因此，我们把这种频率、相位等随时间发生变化的信号称为时变信号。我们通过观察MicroMotion公司生产的型号为CNG050微弯型科氏质量流量传感器的实际工作过程发现：
1）当流量较稳定时，两路信号的相位随时间会发生缓慢而微小的变化，这种变化是随机的没有规律的，但变化的幅度并不像文献[1]中建立的信号模型那样大，一般不超过给定相位的1%。在小流量测量中，如1kg/min的流量对应的相位差为0.020左右，此时相位波动不超过0.00020，变化范围很小，但要提高小流量测量精度这种波动是不能忽略的。
2）当流体密度确定后，传感器信号的频率也会发生变化，但变化范围比相位要小很多，最大不超过流量管振动频率的万分之一。同时信号频率反映流体的密度，因此算法必须适用于测试不同的流体，也需要跟踪信号的频率。
为此，改进时变信号模型如式（1）和（2），更加逼真地模拟实际科氏流量计信号，便于算法仿真。
式（1）和（2）表示信号的相位、频率在变化，每一时刻的值是前一时刻的值加一个随机数，其中eΦ（n）和eω（n）为零均值、正态分布、方差为1且互不相关的白噪声，σΦ和σω分别控制eφ（n）和eω（n）的变化幅度，当信号变化缓慢时，两者减小，当信号突变时，两者增大。λΦ和λΦ分别控制Φ（n）和ω（n）的变化幅度，相位变化幅度应低于给定相位的1%，频率变化低于振动频率的0.01%，这样比较符合实际情况。
3 算法原理及推导
3.1 多抽一滤波
为了增强对噪声的抑制，先用16kHz较高的采样频率对科氏流量计的输出信号进行采样，然后用多抽一滤波器进行抗混叠滤波和抽取。多抽一滤波器分为两级[4]，第一级为2抽1，使实际的采样频率从16kHz降低到8kHz，主要目的是减少数据量。第二级为4抽1，采样频率降低为2kHz。同时采用30阶FIR低通滤波器，不仅保证线性相位，而且在实际的实现中，可以只对抽取的点进行滤波，然后再抽取，这样可以减少计算量节省时间。多抽一滤波器的系数在确定截止频率后通过计算机辅助设计的方法得到。仿真结果表明该方法由于尽可能多地获取了原信号的信息，所以比单纯用2kHz采样、滤波所得的效果要好。
3.2 自适应格型陷波滤波器
自适应陷波滤波器参数可以根据信号特征收敛并可估算信号的频率。采用的格型IIR陷波器[1]是由全极点和全零点两个格型滤波器级联而成，传递函数为：
（3）
为了减少计算负担，通过将零点固定在单位圆上，使得只调整一个参数就能达到自适应陷波的目的。将零点固定在单位圆上，即令k1=1，k0在经过一段时间自适应后收敛到−cosω，ω是信号的归一化频率，α决定陷阱的宽度，k0使用Burg算法[1]进行自适应调整。
由于科氏流量计流体的密度反映为频率的变化，需要及时跟踪流体信号的频率变化。通过大量仿真研究发现，通过调整ρ和λ的终值，适当地加大陷波器陷阱的宽度，便能在保证精度的同时实现对频率变化的跟踪，ρ和λ计算公式如式（4）和（5）所示。
（4）
（5）
格型自适应陷波滤波器的计算量大大降低，且参数少调整方便。调整ρ和λ的终值以及变化的步长就能方便的跟踪频率的变化，同时亦能达到很高的精度。
3.3 SDTFT递推算法及测量相位差原理
3.3.1 SDTFT递推算法
离散时间序列的傅里叶变换（DTFT）为：
（6）
DTFT是从第一个采样点开始通过不断增加计算的序列长度来实现指定频率处傅里叶系数的计算，如果信号在一段时间内恒定不变，这种算法是可行的。但是，无法用于时变信号。时变信号的每个采样点都包含着相位变化的新信息，DTFT将相位变化的新旧信息全部混淆叠加在一起，对相位的变化根本无法灵敏的反映。因此，我们提出滑动DTFT来处理时变信号。
给所观测的信号加一个N点的时间窗，矩形窗是最简单的时间窗，并让这个时间窗随着采样点数的增加不断向前滑动，如图1所示。随着窗函数的滑动，在每个采样点计算N点有限长序列的傅里叶变换即为滑动的或滑动窗的DTFT（SDTFT）。面向时变信号时，计算新采样点的傅里叶系数时仅利用的是当前采样点之前的N点（N是可以改变的），更新新的相位信息并摒除旧的相位信息，这样时间窗随着新采样点不断向前滑动，计算的相位差才能跟踪上实际相位差的变化。
图1 N点滑动时间窗
如图1（a）所示，对于观察信号x（t），设在m时刻采样得到N个采样数据x（0），x（1），…，x（N–1），首次构成N点有限长序列，其离散时间傅里叶变换为：
（7）
式中：ω为数字角频率，单位为rad，t表示采样点的序号。
图1（b）所示在m+1时刻，得到新的采样点x（N），则该点与之前的N–1点重新构成一个N点有限长序列，该序列在处的离散时间傅里叶变换为：
以此递推，当新的采样点与其之前的N–1个采样点组成第k个N点时间窗即采样点序号为（N+k–1）时，该序列在ω处的傅里叶变换如式（9）所示。
式（9）即为SDTFT的递推算法的递推公式。可见，每采入一点新的信号，虽然需计算N点傅里叶变换，但通过递推公式并没有增大计算量，比滑动Goertzel算法计算量小，可以满足科氏流量计实时性要求。同时，每计算一个采样点的傅里叶系数始终是N点叠加的计算结果，不存在序列不断叠加溢出的问题，非常利于实际系统的实现。
3.3.2 窗长度N的选取
科氏流量计传感器信号是正弦信号，具有周期性性质。同时，格型滤波器估计频率精度虽然较高，但仍然与真实频率有偏差。这样，应用DTFT及SDTFT仿真计算周期信号在有偏频率下的傅里叶系数时会出现如图2的现象。
图2 计算有偏频率下的相位差
图2中点划线为真实相位差值，实线为两种方法的估计值。在图2的上图中，可以看出由DTFT计算的每个采样点输出的相位差会发生波动，但波动是偏离真实值的，且真实值偏上的部分比偏下的部分幅度小，这样平均处理后得到的结果会比真实值偏小；而图2中下图所示的SDTFT计算结果中，每个采样点输出的相位差是在真实值上下波动的，虽然比上图中波动的幅度要大，但是真实值上下波动幅度相当，这样平均处理后会比DTFT更接近真实值，精度更高。出现图2仿真结果是因为鉴于处理信号的周期性，我们根据格型滤波器估计的频率值选择SDTFT的窗函数长度N尽量接近信号周期的整数倍，同时SDTFT对整周期的要求要远远低于SDFT[6-7]的要求。因此SDTFT的时间窗函数长度的选择要针对处理信号特点选择。
同时，窗长度N的选取还要根据实际信号的具体变化灵活选取，如果信号的相位变化比较缓慢，可以将N增长，不仅能跟踪上变化，同时点数多可以提高计算精度；如果信号的相位差变化快速，可将N缩短，增加跟踪速度，但不可避免地牺牲了精度，此时可以通过改变窗函数的形状如汉宁窗等来提高计算精度。
3.3.3 SDTFT递推算法计算相位差
由于科氏质量流量计传感器信号为正弦信号，因此可进行计及负频率的修正[2]，减小频谱中负频率成分的影响，增加计算相位差的精度，缩短收敛过程。具体推到公式如下：
根据式（12）得到相位差和时间差后，即可根据仪表系数得到瞬时流量和累积流量。本文测试相位差的方法用于时变信号时不仅能跟踪微小缓慢的变化，而且跟踪速度和精度均优于滑动Goertzal算法。
4 MATLAB仿真结果
型号为CNG050的科氏质量流量传感器，满管振动时传感器信号基频为188.64Hz，因此仿真时信号频率采用188.64Hz，着重仿真小流量对应的相位差。
根据时变信号模型产生的相位随采样点数变化的曲线以及生成的时变信号波形如图3所示。
图（3）、图（4）仿真参数为：
（a）相位变化
（b）时变信号波形
图3 按照时变信号模型产生的时变信号
图4 本文方法跟踪相位差变化效果
图4所示即相位在[0.009940，0.010140]内变化时，本文算法的跟踪效果图，可见波动幅度为0.01的1%时，本文方法仍具有很好的跟踪速度和跟踪精度，而SGA算法已经无法跟踪此时相位的微小变化了。
图5所示的是相位在[0.0080，0.0150]内变化时，本文算法与SGA算法跟踪相位变化的效果对比，可以看出，相位变化超过0.010的70%时，SGA才能跟踪上变化，但本文算法的跟踪速度和精度都优于SGA。
图5 SDTFT与SGA跟踪相位比较
[0.010，0.20]内几个相位的仿真结果数据如表1所示，仿真参数同图6。从表1中可以看出，本文方法具较高的精度。
5 系统研制
本课题组研制了相应的数字式科氏质量流量变送器系统，一方面为激振器提供激振信号，另一方面对两路传感器信号进行处理，得到流量信息。
5.1 系统硬件组成
系统硬件基于TI公司TMS320F28335浮点DSP，由驱动模块、信号调理采集系模块、脉冲和4～20mA电流输出及人机接口组成。
5.2 系统软件研制
系统软件设计采取模块化设计方案，将完成特定功能或者类似功能的子程序组合成功能模块，主要功能模块如图6软件总体框图所示。
图6 系统软件总体框图
5.2.1 主要模块
初始化模块负责系统内可编程器件和全局变量等的初始化，包括TMS320F28335初始化、全局变量初始化、算法初始化和LCD初始化等。F28335的初始化[8]主要包括系统配置、用到的28335内部集成功能模块的初始化，如看门狗、多通道缓冲串口McBSP、DMA、系统外部接口XINTF和数字复用口等以及外部可编程器件如CodecCS4271芯片。全局变量以及算法初始化主要包括反映系统工作状态的状态变量和算法计算的初始变量等。LCD初始化则包括控制引脚初始化、写命令字和初始化显示内容。
中断模块由外设中断扩展控制器来处理所有片内外设和外部引脚中断的优先级以及中断的响应。本系统采用了三个中断，DMACH1接收中断、DMACH1发送中断、定时器0的周期中断触发。Codec启动后开始采传感器的信号并进行AD转换，之后送到DSP的McBSP的接收引脚上，通过DMA传输到内部RAM中，当RAM放满数据后产生DMACH1接收中断，中断响应后将数据放于外部RAM一个长的循环数组中供算法调用，并修改DMA的目的地址，为下次传输做准备；定时器0的周期中断触发用于定时刷新LCD的显示值。
算法模块内包含了系统所进行的大部分数学运算子程序，主要有信号处理算法、瞬时流量计算和累计流量计算以及温度补偿算法。信号处理算法在主监控程序中调用，先调用多抽一算法程序进行抽取，随后调用格型滤波估算频率，之后再利用负频率修正的SDTFT递推算法计算信号的相位差，平均处理后，再结合仪表系数计算当前的瞬时流量。还需根据特定公式对结果进行温度补偿。测量结果包括瞬时质量流量、累积质量流量、密度和温度，由LCD显示。
主监控程序[9]是整个系统的总调度程序，实现仪表所要求的功能。系统上电复位后，立即进行初始化，完成对系统及各模块的初始化工作，然后监控程序开始依次查询标志寄存器的有关位，如采样数据标志位、频率跟踪标志位、流量计算标志位、通信标志位等，如果有标志位被置位，则调用相应的子程序，并清除标志位为下一次操作做好准备。一次查询完后，监控程序进入等待状态，直到有关中断发生时开始进入下一次查询，周期循环，实现仪表流量的实时检测。
5.2.2 关键实现技术
1）TMS320F28335是TI公司推出的一款新的浮点运算DSP芯片，由于负频率修正的SDTFT算法计算小相位差的关键在于保证公式（20）中分子得到的数据的有效位数，经过研究发现用32位float型变量类型在小流量时只能保证3位有效数字，使得计算精度降低，因此我们使用64位longdouble型变量类型；
2）本文提出的一套数字信号处理算法有大量的变量，且变量采用64位的数据类型，这就需要较大的存储空间，因此外扩了64k的SARAM，为了提高计算速率，在存储器映射设置时将全局变量全部映射到外部存储器中，将使用频率高的局部变量和程序映射到内部RAM中，使得算法成功实现；
3）经过测试F28335定点CPU计算外部RAM中一点数据的两级多抽一算法时间为116μs，计算格型滤波算法的时间为83μs，计算SDTFT递推算法时间为221μs，再加上对结果的补偿及平均处理，得到一点流量的时间不超过500μs，因此CodecAD采样频率设置为16kHz，多抽一后为2kHz，在保证精度的情况下完全实现了连续采样数据的实时计算；
4）AD设置较高的采样频率，这就要求传送大量的数据，一点一点传送显然不满足实时性的要求，为此我们利用了TMS320F28335的多通道缓冲串口McBSP与CodecAD通讯，运用DMA传输大量的数据而不降低CPU的使用效率；
5）在系统软件的初始化模块中对各个部分的初始化顺序有严格的要求，否则系统无法正常运行或者不稳定。例如在F28335系统配置之后，要先初始化通用输入输出口；又例如要先初始化Codec芯片再初始化McBSP和DMA模块等。
6 系统测试
考虑到在实际流量标定中，很难产生变化规律已知的时变信号，所以，我们采用以下方式对本文提出的时变信号处理方法进行测试。首先根据时变信号模型，利用MATLAB产生已知参数的时变信号数据；其次，将信号数据导入Fluke282任意波形信号发生器，由信号发生器产生两路存在相位差且时变的信号；再将信号发生器的两路输出分别接至信号处理系统的两个输入通道，运行DSP程序，对信号进行处理。相位差在[0.010，20]范围、变化为1%的测试结果如表2所示。
由表2可以看出小相位差的相对误差大于0.1%，结果并没有仿真时精度高，其误差来源为：1）科氏流量计变送器系统采用Codec芯片采集数据，其AD位数为13位，从而限制了计算精度；2）MATLAB产生的信号是64位、双精度浮点数，但Fluke282信号发生器的数据是12位，这造成了截断误差。
7 结论
1）面向时变信号，提出由多抽一滤波、自适应格型陷波滤波、负频率修正的SDTFT递推算法组合而成的一套数字信号处理算法，由仿真结果可以看出该算法能在每个采样点上输出傅里叶系数，实时性好，能跟踪信号微小的变化，在小流量时仍具有较高的精度，性能优越，且算法计算量小，不存在数值溢出；
2）将整套算法在TMS320F28335DSP搭建的数字式科氏质量流量变送器系统上实时实现，得到了较好的效果。可见这种处理方法能够实时获得信号间的相位差和时间差，可以提高科氏流量计的动态响应速度，满足一些测量场合的需要，具有较强的实用性。
参考文献
[1]徐科军，倪伟，陈智渊.基于时变信号模型和格型陷波器的科里奥利质量流量计信号处理方法[J].仪器仪表学报，2006，27（6）：596-601.
[2]李祥刚，徐科军.科氏质量流量管非线性幅值控制方法研究[J].电子测量与仪器学报，2009，6：85-99.
[3]张海涛，涂亚庆.计及负频率影响的科里奥利质量流量计信号处理方法[J].仪器仪表学报，2007，3（3）：540-541.ZHANGHT，
[4]TUYQ，ZHANGHT.MethodforCMFsignalprocess-ingbasedontherecursiveDTFTalgorithmwithnegativefrequencycontribution[J].IEEETransactionsOnInstru-mentationAndMeasurement，2008，57（11）：2647-2653.
[5]DERBYHV，BOSET，RAJANS.Methodandapparatusforadaptivelineenhancementincoriolismassflowmetermeasurement[P].USPatentNo.5555190，Sep.10，1996
[6]JACOBSENE，LYONSR.TheslidingDFT[J].IeeeSig-nalProcessingMagazine，vol.20，pp.74-80，Mar2003.
[7]JACOBSENE，LYONSR.AnupdatetotheslidingDFTieeesignalprocessingmagazine，2004.
[8]TexasInstruments.TMS320F28335，TMS320F28334，TMS320F28332，TMS320F28235，DigitalSignalCon-trollers（DSCs）DataManual.February，2008.
[9]徐科军，于翠欣，苏建徽，等.基于DSP的科式质量流量计信号处理系统[J].仪器仪表学报，2002，23（2）：
客服1