基于MATLAB GUI编程的血流/温度信号时频分析软件开发

唐元梁，贺缨

大连理工大学 能源与动力学院，辽宁 大连 116024

引言

微循环是人体血液循环系统的重要组成部分，其基本功能是进行血液和组织液的物质交换，因此健全的微循环功能是人体重要脏器功能正常运行的重要保障。许多疾病的发生，包括高血压、糖尿病及心脑血管疾病等，均与微循环障碍有关。若能及时检测出微循环改变的危险信号，将对微循环疾病的早期诊断和治疗有极大帮助。

通过对人体皮肤血流信号进行频谱分析，研究者发现其包含多个波动特征频率[1-4]（图1）。虽然这些波动并非是严格周期性的，但很显然也不是随机出现，表现出类周期性的特点。通过干预交感神经活动、内皮活动等，研究者发现这些波动特征频率与微循环的调节活动有关，这些调节活动及其活动周期范围如下：内皮调节（50～105 s）；神经调节（20～50 s）；肌源性调节（7～20 s）；呼吸（2～7 s），心跳（0.5～2 s）。因此，若能通过滤波提取出与各调节活动相对应的波动分量并进行数据分析，可以更加准确地判断微循环调节功能病变的发生原因。

图1 人体皮肤血流信号的波动特征

由于皮肤温度受皮下微循环血液灌注的影响，因此能够间接地反映皮下微循环的状况，且相较于血流信号的采集（一般采用激光多普勒血流仪等），体表温度采集更容易、测量稳定性更好。而且实验分析显示[5-7]，人手指尖等部位的皮肤温度也表现出与血流信号类似的波动特征，且两者具有一定的波动相关性。因此，综合考虑微循环调节活动的类周期性特征、皮肤温度与血流的波动相关性以及热信号采集方面的优势，Frick等[8]指出，皮肤温度的波动特征可以用来标记微循环的调节功能。近年来微循环功能的热波动检测逐渐引起的研究者的关注。Smirnova等[9]发现在单侧手接受冷水刺激时，健康人群的对侧手皮肤温度在内皮、神经、肌源性调节频段的波动幅值均减小，冷刺激结束后波动幅值迅速回升至正常；而二型糖尿病、糖耐量受损人群在结束冷刺激后，3个频段内的温度波动幅值均没有恢复正常，表明二型糖尿病组、糖耐量受损组的皮肤微循环调节功能均出现了障碍。研究者发现指尖局部热刺激时健康人群的皮肤温度在肌源、神经以及内皮调节频段的波动幅值都有明显的增加[10-12]，这反映了热刺激诱导的血管舒张调节变化；而二型糖尿病患者的指尖在热刺激前后其内皮调节频段的温度波动幅值并没有显著变化，表明后者可能存在内皮调节功能障碍。

微循环功能热波动分析方法可以应用于相关软、硬件产品的开发。FM Diagnostics公司与俄罗斯科学院连续介质研究所共同研发了Microtest仪器[12]，见图2，该仪器主要包括了机身、探头两部分。探头部分尺寸约20 mm×30 mm×10 mm，温度传感器（HRTS-5760，Honeywell International，Inc.，Morristown，NJ，USA） 安装在探头圆孔（直径约8 mm）的中央，传感器外周有4个小型电阻加热器，探头内填充低热导率材料来减少环境变换的影响。经信号放大和A/D转换后温度测量精度可达0.001℃。将人体待检测部位皮肤紧贴Microtest温度传感器，启动仪器后可以对皮肤进行变功率加热和皮肤温度采集，然后通过配套的软件自动对获取的皮肤温度信号进行时频分析，并根据不同功率热刺激下皮肤温度波动幅值的变化情况对皮下微循环功能进行评价。Sagaidachnyi等[7,13-15]开发了TBF Converter软件，该软件提供了信号的Fourier分析和小波分析、PPG信号包迹线提取功能，并能根据不同模型将血流和温度信号进行相互转化（图3）。

图2 Microtest仪器

利用血流或温度信号的波动特征对微循环调节功能进行评价时，均首先需要对信号进行频谱分析。信号分析中，傅里叶变换是常用的方法，但经过傅里叶变换后时域信息丢失，因此傅里叶变换只在频域里有局部分析能力，采用标准傅里叶变换不能得到信号波动的时变特征。短时傅里叶变化和小波变换均能保留时域信息，但相较于窗口宽度固定、只有单一的分辨率的短时傅里叶变化，小波变换采用可变的时间窗，即低频信号可采用较宽的时间窗，而高频信号采用较窄的时间窗，因此可以在时、频两域都具有表征信号局部特征的能力。小波变换非常适合于探测正常信号中夹带的瞬态反常现象并展示其成分，誉为信号分析的“显微镜”。作为一种具有多分辨率特点的时间—尺度分析方法，小波变换是分析微循环信号（或皮肤温度信号）不同频率波动成分及其变化特征的非常有用的数学工具。

从目前情况来看，尽管基于热波动分析方法对皮下微循环调节功能进行检测的前景非常可观，但相关的软、硬件的发展还处于比较早期的阶段，距临床应用还有较大距离。目前已有的相关产品数量较少，功能较为简单。如Microtest仪器对于微循环检测的模式只有加热刺激这一种，而且其配套软件只能处理该仪器获取的数据，不能够从外部导入；TBF Converter软件侧重于信号分析和转换，没有信息挖掘功能，不能直接用于微循环功能的评价。尽管小波分析方法非常适合于微循环信号分析，但其实现过程比传统的傅立叶分析更复杂一些，如小波函数的选取、尺度与频率的转换以及边缘效应等问题均可能会对研究者带来困扰。虽然MATLAB平台提供的Wavelet Toolbox工具箱能极大地简化用户进行信号小波分析的过程，但该工具箱仅包含基本的小波变换功能，后续的数据处理，如不同频段上信号波动分量的提取、波动幅值的分段统计，不同数据间的波动相关性分析，以及数据处理结果的图像/表格显示等仍需用户自己编程实现。

图3 TBF Converter软件界面

因此，为便于使用小波分析方法对微循环血流/温度信号进行波动分析，本文基于MATLAB平台的guide工具箱开发了一个有图像用户界面的应用软件，期望对血流/温度信号的小波分析、微循环功能的检测提供帮助。

1 生物信号的小波分析方法

1.1 一维连续小波变换

小波分析将信号分解成一系列小波函数的叠加，而这些小波函数都是由一个母小波函数经过平移与伸缩得来的（图4）。母小波函数φ(t)经伸缩、平移后可得小波序列φa,b(t)：

其中a为尺度因子，b为平移因子，两者均为实数。对信号x(t)进行连续小波变换后得到小波系数Wx(a,b)：

式中T1、T2分别为时域信号x(t)的起、止时间，Wx(a,b)为信号x(t)的小波系数，φ∗为φ的共轭函数。尺度因子a与频率f呈倒数关系：

其中fcen为母小波函数φ(t)的中心频率，Tsam为信号x(t)的采样周期。在信号分析领域，Morlet小波是最常用的小波函数，其具有良好的时-频定位性质，时域下Morlet复小波函数的简化形式为：

图 4 一维连续小波变换图示

1.2 信号频谱特征分析

类似于傅立叶分析，小波分析也可以提取信号的频谱特征。小波系数的模方表征了信号的小波功率谱，将小波功率谱进行时间平均可得到某时间段上信号的小波全谱M(a)：

小波全谱M(a)反映了在尺度a上信号的波动强度，根据频率与尺度a的关系[式(3)]，小波全谱图的局部峰值代表了信号的波动特征频率（或周期）。

1.3 波动相关性分析

信号x(t)、y(t)的整体互相关性可以利用如下的互相关函数计算：

1.4 波动分量提取

利用信号的小波系数Wx(a,b)进行小波反变换，可以提取出原始信号在某一频段上的波动分量，如对频段[f1, f2]上信号x(t)的波动分量x1(t)进行提取：

其中[a2, a1]为频段[f1, f2]对应的尺度因子范围。提取出波动分量x1(t)后，可以利用其标准差来表征该频段上信号波动幅值的大小。

1.5 信号时频分析软件的开发

为了便于对血流/温度信号进行基于小波变换的时频分析，本文通过MATLAB GUI编写了可视化的数学程序，并生成可执行的软件（exe文件）。本文所设计开发的软件主要包括了文件模块、小波分析模块、检测结果模块以及相关性分析模块4部分，见图5。

图5 软件模块

文件模块主要功能是导入血流、温度等信号并进行显示；相关性分析模块主要功能是计算导入的血流、温度信号的互相关性；小波分析模块可以对导入的信号进行小波变换，并提取出不同频段的信号波动分量幅值变化信息；若导入的血流/温度信号是血管反应性测试实验中获取的数据，如反应性充血、热刺激和冷刺激等，则根据信号波动分量的幅值变化信息，可以在检测结果模块中统计出信号在不同实验阶段的波动幅值变化情况。

在MATLAB软件的命令行中输入“guide”，便可创建新的GUI。将编写的所有GUI添加到一个工程目录下并编译通过后，便可以打包生成可执行程序，双击该程序便可运行。不过需注意的是MATLAB GUI编程所得的exe软件的正常运行依赖于MATLAB数学函数库，因此若计算机上未安装MATLAB软件，需要安装MCR（MATLAB Compiler Runtime）编译器。

2 软件使用介绍

2.1 信号小波变换和波动幅值提取

软件的小波分析模块界面如图6所示。其中区域①显示了导入的血流时间序列信号，该信号为利用ADI公司的ML191激光多普勒血流仪采集，采样周期0.001 s；信号采样周期默认为1.0 s，小波分析尺度默认为1～100，根据信号的实际采样周期和小波分析尺度要求，在“采样周期”和“分析尺度”后方的文本编辑框内可修正相应的数值；点击“小波时频分析”按钮，软件对导入的温度信号进行小波分析，并对分析结果进行显示；区域②可以显示某段时间内信号的小波系数图和小波系数模图，图像模式可在左上方的弹出菜单选择，时间参数可在“开始时间”和“分析时长”后方的文本编辑框内设置，图中显示的是区域①所示信号在15～21 min内（对应区域①中两条蓝线间的信号部分）的小波系数图；区域③显示了信号的小波方差图；区域④的3条曲线由上到下分别为所输入信号在内皮、神经和肌源性频段的波动幅值变化，默认的尺度因子变化范围与人体微循环3个主动调节活动的频段相一致，如需调整则可分别在“最大尺度”和“最小尺度”方的文本编辑框内修改。

2.2 检测结果的分析

在某些机械、温度刺激条件下，正常皮肤微循环的调节活动强度会发生变化，而有微循环功能障碍的人群其调节活动响应可能会出现异常。这些调节反应的状况可以通过皮肤血流/温度信号波动幅值的变化进行评估。

本软件在“检测结果”模块中，在对不同刺激实验（热刺激、对侧手冷水刺激、反应性充血）获取的温度信号进行小波分析后，将提取的波动幅值进行了分段统计；根据不同实验阶段的温度波动幅值变化并参考已有文献中的温度波动变化指标，可以对（对侧手冷水/热刺激下）皮下微循环状态进行定性评价。软件的检测结果模块的界面，见图7。图中显示的是一人手指尖热刺激实验获取的皮肤温度信号在不同实验阶段的平均波动幅值变化情况。该结果显示Stage 1阶段各频段的平均温度波动幅值均高于Stage 2阶段，表明微循环调节活动强度变化对热刺激敏感，调节功能正常。

除温度刺激实验检测结果分析外，该模块还提供了反应性充血实验数据的分析功能。利用大鼠反应性充血实验，本研究小组发现二型糖尿病组大鼠在血管闭塞阶段内皮、神经调节频段的皮肤温度波动幅值小于对照组大鼠且差异明显[16]；刺激结束后，二型糖尿病组大鼠神经、肌源性频段的皮肤温度波动幅值大于对照组大鼠且差异明显。该实验结果对于评价人体糖尿病微循环障碍有一定价值，但还需大量的人体实验进行验证。

图6 小波分析模块

图7 检测结果模块

2.3 信号间的相关性分析

图8 互相关分析模块

软件的相关性分析模块的界面，见图8。其中血流、温度图的横坐标为时间，单位s；小波互相关谱图的横坐标为波动周期，单位s。图中红色、绿色曲线分别为反应性充血实验获取的相同位置指尖皮肤血流和温度信号，其中皮肤血流信号采用激光多普勒血流仪采集，温度信号采用Agilent温度记录仪采集，精度为0.001℃。经互相关分析后求得两信号的整体互相关系数为0.7385；黑色曲线为两信号在1～100尺度范围内的小波互相关谱，从整体上看两信号在较大尺度上波动相关性较高，表明该温度信号能够表征血流信号的低频脉动特征。

3 总结

利用MATLAB GUI编程，本文基于小波分析方法开发了一款用于微循环血流/温度信号时频分析的软件。在小波分析模块中，软件提供了小波变换功能，可以对信号的频谱特征进行分析；除了基本的小波变换功能外，本软件还能够提取信号在不同频段的波动幅值变化信息，并能够对导入的血流/温度信号进行波动相关性分析。由于在微循环调节活动影响下人体皮肤血流信号呈现出类周期性特征，因此皮肤血流信号波动幅值变化信息能够反映皮下微循环调节活动是否出现异常。而当血流和温度信号有较高波动相关性时，皮肤温度信号波动幅值的变化也能用于微循环功能的检测。因此，结合该软件各模块的不同功能，可以为微循环的热波动检测提供较为直观、可信的数据信息。

该软件还需进一步完善。目前我们在检测结果模块中只对信号波动幅值变化进行了分段统计，还需更多量化指标以全面评价微循环。我们将进一步通过动物和临床实验获取更多数据信息，并通过数据分析来逐步优化软件，以适应将来微血管功能评价的临床应用需求，为糖尿病足等的早期诊断提供帮助。

[参考文献]

[1] Kvernmo HD,Stefanovska A,Bracic M,et al.Spectral analysis of the laser Doppler perfusion signal in human skin before and after exercise[J].Microvasc Res,1998,56:173-182.

[2] Bracic M,Stefanovska A.Wavelet-based analysis of human blood- flow dynamics[J].B Math Biol,1998,60(5):919-935.

[3] Kvernmo HD,Stefanovska A,Kirkebøen KA,et al.Oscillations in the human cutaneous blood perfusion signal modified by endothelium-dependent and endothelium-independent vasodilators[J].Microvasc Res,1999,57(3):298-309.

[4] Soderstrom T,Stefanovska A,Veber M,et al.Involvement of sympathetic nerve activity in skin blood flow oscillations in humans[J].Am J Physiol Heart C,2003,284(5):1638-1646.

[5] Frick P.Wavelet-based correlations of skin temperature and blood flow oscillations[J].Cardiovasc Eng,2008,8(3):185-189.

[6] Bandrivskyy A,Bernjak A,Mcclintock P.Wavelet phase coherence analysis: Application to skin temperature and blood flow[J].Cardiovasc Engin,2004,4:89-93.

[7] Sagaidachnyi AA,Skripal AV,Fomin AV,et al.Determination of the amplitude and phase relationships between oscillations in skin temperature and photoplethysmography-measured blood flow in fingertips[J].Physiol Meas,2014,35(2):153-166.

[8] Frick P,Mizeva I,Podtaev S.Skin temperature variations as a tracer of microvessel tone skin temperature variations as a tracer of microvessel tone[J].Biomed Signal Proces,2015,21:1-7.

[9] Smirnova E,Podtaev S,Mizeva I,et al.Assessment of endothelial dysfunction in patients with impaired glucose tolerance during a cold pressor test[J].Diab Vasc Dis Res,2013,10(6):489-497.

[10] Parshakov A,Zubareva N,Podtaev S,et al.Local heating test for detection of microcirculation abnormalities in patients with diabetes-related foot complications[J].Adv Skin Wound Care,2017,30(4):158.

[11] Podtaev S,Stepanov R,Smirnova E,et al.Wavelet-analysis of skin temperature oscillations during local heating for revealing endothelial dysfunction[J].Microvasc Res,2015,97:109-114.

[12] Parshakov A,Zubareva N,Podtaev S,et al.Detection of endothelial dysfunction using skin temperature oscillations analysis during local heating in patients with peripheral arterial disease[J].Microcirculation,2016,23:406-415.

[13] Sagaidachnyi AA,Usanov DA,Skripal AV,et al.Correlation of skin temperature and blood flow oscillations[J].Proceedings of SPIE,2012,8337(4):9.

[14] Sagaidachnyi A,Usanov DA,Skripal AV,et al.Restoration of finger blood flow oscillations by means of thermal imaging[A].11thInternational Conference on Quantitative Infra Red Thermography[C].2012.

[15] Sagaidachnyi AA,Usanov DA,Skripal AV,et al.Thermographybased blood perfusion imaging in hands : spectral amplification and time shift[A].Quantitative Infrared Thermography[C].2014.

[16] 唐元梁,糖尿病微血管病变热检测方法的机理研究[D].大连:大连理工大学,2017.