基于LabView和Matlab混合编程的脑MIT数字鉴相方法的实现

【摘要】  目的： 为实现用于脑磁感应成像系统中的数字鉴相方法，寻求一种高精度的相位检测手段. 方法： 本研究采用LabView对数据采集卡编程，控制激励信号的产生和检测信号的采集；采用Matlab进行快速傅立叶变换法(FFT)法，相关法和经典法的编程，用于实现数字鉴相. 结果： 得到了步进0.1°条件下的数字鉴相结果，构建了一整套包含信号产生、采集和处理的鉴相系统，最终，实现了FFT法和相关法0.1°的相位差分辨率. 结论： 数字鉴相方法能够应用于脑磁感应断层成像技术（MIT）检测系统中，对高精度的系统测量提供了理论依据.

【关键词】  脑磁感应成像；相位检测；导电性；傅里叶分析；相关法

　　在脑磁感应断层成像技术(magnetic induction tomography, MIT)中，激励信号与检测信号之间的相位差产生源于磁场中脑组织的存在，同时，对组织电导率的求取可以转化为相位差的求取［1-3］. 因此寻求一种高精度的相位检测方法对系统设计具有重要意义. Watson等［4］和Scharfatter等［5］多采用模拟鉴相器的方法进行MIT方面的研究，而采用软件数字鉴相方面仍处于起步阶段. 本研究针对模拟鉴相器容易引入噪声，精度难以提高的现象，拟采用数字鉴相方法对MIT监测系统进行改进，旨在实现一种低误差、高精度的软件鉴相方法.

　　1材料和方法

　　1.1材料

　　本系统采用了基于LabView的虚拟仪器接口编程方法［6］实现信号产生和采集. 数据采集卡使用美国国家仪器公司(National InstrumentsTM, NI)的NI?6115，最高采样频率10 MHz，12位精度，输出电压范围，最大双通道同步数模转换器(digital and analog convertor, DAC)的输出频率为2.5 Msamps/s，程序采用7.0版本编写. 数字鉴相算法采用美国Mathwork公司推出的Matrix Laboratory(Matlab)实现.

　　1.2方法

　　1.2.1系统设计

　　MIT检测系统中检测信号与激励信号间存在由于生物组织所产生的相位差，它们的模拟利用数据采集卡编程实现，对相位的三种检测在PC内部完成，以便实现方便控制的检测方法. 具体方法如图1：

　　1.2.2信号产生

　　MIT系统中，初级磁场的强弱对信号检测有着非常大的影响，过大的初级磁场会使检测信号被淹没，因此考虑使用低频信号进行处理. 信号产生由NI?6115通过LabView编程实现，同步发生两路200 KHz正弦信号（幅度、相位、波形、偏移量可调），由于NI?6115双通道DAC的最高输出频率为2.5 Msamp/s，因此每周期最多只能产生12个点. 信号的产生采用LabView中级模拟输出函数. 包括通道设置模出操作，启动带缓冲的模出操作及清除模出操作. 程序运行时，设置模出操作会产生一个任务编码(taskID)和出错信息簇，其余的模拟输出接受taskID以识别操作的设备和通道，并且在操作完成后输出一个taskID，该参数形成了采集函数之间的一个关联数据.

　　将产生的信号进行谱估计，发现产生的信号噪声成分很多，通过，信号频率能量占总能量的59.48%，因此考虑对信号进行滤波处理. 能量计算公式如下：ε=A2〖〗∑N〖〗i=0X2i (1)
(1)式中的A代表有用信号的幅度，N代表采样点数，Xi代表频谱中各个频率成分信号的幅度.

　　1.2.3信号采集

　　数据采集由另外一台PC完成，通过两根双同轴电缆头连接线(bayonet nut connector, BNC)传输线连接端口，采样速率均为10M/S（NI?6115最大的采样频率），每周期采集50个点，共采集40周期数据进行存储. 进程包括通道设置模入操作，启动带缓冲的模入操作，从被分配的缓冲读取数据，以及清除模入操作在计算机中分配的缓冲、释放所有数据采集卡的资源. 其中还包括一子程序，在其内部可以进行数据处理，然后输出，此处没有对数据进行处理，完成的功能仅仅是判断数据列是否出错，如果出错则退出循环，否则正常进行程序流程.

　　1.2.4数据处理

　　数字处理过程根据不同算法的要求，对数据进行不同处理，快速傅立叶变换法(fast fourier transform, FFT)和相关法采用200阶有限长冲击响应(finite impulse response, FIR)数字带通滤波器，得到待测信号后进行FFT和相关处理；经典法先将数据通过数字带通滤波器，然后将两路信号相乘，提取直流分量，最后计算相位差，绘制曲线. 数字鉴相算法由Matlab编程实现，能够实现采集信号处理和Matlab内部仿真信号处理的功能.

　　信号处理软件采用了三种数字鉴相方法. FFT法中采用了以下结论［7］：X(l)=NA〖〗2ejφ(2)
(2)式中的N代表采样点数，A代表信号幅度，l代表频谱中标志被测信号频率的谱线序列号. 对复数X(l)求反正切即可得到初始相位φ. 两路信号同时进行以上处理并相减即可得到相位差. 具体核心程序如下：
X=fft(x);Y=fft(y);//两路信号进行FFT变换
theta1=angle(X(41))*180/pi//41代表频谱中第41个数据是信号的谱线
theta2=angle(Y(41))*180/pi//对两路信号进行寻谱并求取初始相位
theta=theta2?theta1;//求取相位差
相关法利用相关系数和信号相位差的关系进行编程［8］，结论如下：
Rxy(0)=AB〖〗2cos(φ1-φ0)(3)
φ1-φ0=arccos(2Rxy(0)〖〗AB)(4)
其中A=2Rxx(0)，B=2Ryy(0). 核心程序如下：
rx=sum(x.*x)/N;
ry=sum(y.*y)/N;//求取两路信号的自相关系数
rxy=sum(x.*y)/N;//求取两路信号的互相关系数
sqr=sqrt(rx*ry);
phase2=acos(rxy/sqr)/pi*180;//依照结论求取相位差
经典法较为常用，其基本思想是将两路具有相位差的同频信号相乘，变为一个倍频信号和一个只具有相位差信息的直流分量的叠加，经过低通滤波器以后即可得到相位差信息. 此方法同硬件鉴相器的设计思路一致，属于硬件鉴相的软件实现，但由于直流噪声的存在，使得这种方法的检测结果和理论分析相差较大.

　　2结果

　　2.1预设相位差步进0.1°时的测量结果双路信号发生器由NI?6115实现，产生具有相位差的双路信号，相位差区间为0.1°到0.9°，每次调整信号发生函数，使初始相位差步进0.1°，采集后对数据分别进行数字鉴相，得到如图2所示的检测结果.

　　从图2的结果中可以看出，FFT算法和相关算法能够保持很好的线性度，能较好的满足实验的要求，虽然存在一定程度的误差，但是较好的线性度说明误差值在总体上能够保持恒定，同时，此系统误差由FIR数字滤波器的线性相移带来，因此在结果处理的时候能够比较方便的解决误差带来的影响，便于校正. 经典法检测效果比较差的原因在于直流分量难于提取，并且直流噪声也不容易去除，可以采用两路信号具有一定频差的方法，对整个系统进一步改进.

　　2.2鉴相方法的参数分析为了对结果进行定量的分析，区分方法的优劣，本研究提供了误差水平和非线性度的评估方法. 具体计算方法如下：ε=20log1〖〗Nm∑i=b〖〗i=a|φl-φr| (5)
NL=20log1〖〗(b-a)m∑i=b〖〗i=a|(xi+1-xi-m)|(6)

　　其中ε为误差，NL为非线性度，N为采样点数，φl为理论相位差，φr为实测相位差，m分别为0.1；a, b为预设定相位差的最小值和最大值（对应曲线横坐标的起始和终止值）. 最终，遵循误差越小，非线性度越小，方法越好的原则，对各种方法进行评估. 通过计算得知FFT法，相关法和经典法的误差水平分别为：-24.63 dB，-18.49 dB和14.16 dB；非线性度分别为：-20.4 dB，-21.8 dB，-2.14 dB. 可以看出，FFT法和相关法具有较好的误差水平和线性度. 对于相位检测而言，线性度是衡量的重要标准，由于经典法具有差的线性度，所以在系统的设计中应该不加采用.

　　3讨论

　　本研究通过Labview和Matlab混合编程实现的软件数字鉴相系统有效地解决了测量信号的采集和处理问题，在MIT系统要求的条件下，可以在极短的时间内计算出双路信号的相位差，为成像实验提供依据. 从评价结果不难看出，FFT法和相关法具有更高的线性度，可以应用于MIT检测系统中，同时由于通过程序实现，在今后的改进中便于修改.

　　在数据产生和采集过程中，均发挥了NI?6115的最大效能，实现了双路信号的同步产生和采集，但由于对信号频率的要求较高(200 KHz),所以每周期最多只能产生12个数据点，同时对单周期内采集的点数也产生了限制. 比较Rosell等［2］的系统设计，随同步采集的通道数减少，但实现了更高频率的处理. 鉴于软件编程的可修改性强，对通道数的提高可以在后续工作中完成.

　　数字滤波器的设计过程中，采用线性相位滤波器，这样可以解决滤波器带来的相移问题. 李世俊等［9］采用硬件滤波器，对于其设计是一个难点，具体参数的确定不容易达到理想效果，尤其是高阶滤波器的设计在硬件中基本不可能实现，本研究所采用的数字滤波器阶数达到200阶，不仅使噪声衰减更快，而且带来的相移固定，这些都是硬件滤波器难以实现的.
在模拟鉴相的过程中，存在的诸多噪声很难得到较好的抑制［10］，数字鉴相可通过数字滤波器的选择较好地滤除工频和倍频干扰. 当测量次数趋于无穷大时，随机误差的代数和为零. 因此，由随机噪声产生的误差可采用多次测量求平均值的方法予以消除. 而采用数字的方法可以较容易的得到高信噪比的检测结果.

　　若能够提高数据采集卡的位数和采样频率，同时降低量化误差，进一步采用相干平均等噪声抑制方法对软件进行优化，将有望实现0.01°或更高的精度，使其在脑MIT领域中发挥更大的作用.

【】
  　　［1］王聪,秦明新,董秀珍,等. 磁感应方式电导率测量基础研究［J］.医学物杂志, 2004,21(3)：182-185.

　　［2］Rosell J, Merwa R, Brunner P, et al. A multifrequency magnetic induction tomography system using planar gradiometers: Data collection and calibration［J］. Physiol Meas, 2006,27:271-281.

　　［3］Griffiths H. Magnetic induction tomography Measurement［J］. Sci Technol, 2001, 12: 1126-1131.

　　［4］Watson S, Williams RJ, Griffiths H. Frequency downconversion and phase noise in MIT［J］. Physiol Meas, 2002, 23:189-194.

　　［5］Scharfetter H, Casanas R, Rosell J. Biological tissue characterization by magnetic induction spectroscopy (MIS): Requirements and limitations［J］. IEEE Trans Biomed Eng, 2003, 50: 870-880.

　　［6］清华大学电机系虚拟仪器实验室,北京中科泛华测控技术有限公司编译LabView用户指南(LabView User Manual January 1998 Edition)［M］. 北京：2000,28:1-7.

　　［7］郑南宁. 数字信号处理［M］.西安:西安大学出版社,1997:133-179.

　　［8］曹宇亚. 介质损耗带点检测数字化处理方法的研究［J］.高压电器,2000,3:17-19.

　　［9］李世俊,秦明新,董秀珍,等. 非接触磁感应脑阻抗断层成像系统设计［J］.中国医学物理学杂志,2003,20(1):55-58.

　　［10］秦明新,李世俊,董秀珍,等. 非接触磁感应脑电导率断层成像系统实验系统研究［A］//中国科协2003年学术年会文集［C］. 北京:中国技术出版社,2003:314.