本文介绍了一种界面友好可与医院联网的便携式强迫振荡呼吸阻抗测试系统。该系统采用STM32产生单频率或复频率的正弦振荡信号，经放大电路放大后驱动扬声器产生正弦振荡气流于受试者口腔中，应用STM32的模数据转换器（ADC）模块分别测量由压力传感器和流量传感器得到的模拟电压信号，再经运算后将得到的参数显示在触摸显示屏上，同时可通过GPRS将数据上传到上位机进行保存和显示。最后使用模拟肺和征集志愿者等方式来评价该系统的可靠性。测试结果表明该系统稳定可靠，可以达到监测呼吸阻抗状况的目的。
引用本文： 周垂柳, 王武, 谢联昇, 曾碧新. 一种便携式强迫振荡呼吸阻抗测试仪的设计与实现. 生物医学工程学杂志, 2016, 33(5): 951-957. doi: 10.7507/1001-5515.20160153 复制
引言
慢性阻塞性肺疾病(chronic obstructive pulmonary disease，COPD)和哮喘(asthma)等是常见的呼吸系统疾病。COPD是一种慢性肺部疾病，主要表现为呼出气流受限不完全可逆和呼吸困难[1]。在我国，COPD在15岁以上的人群患病率为3%，全国患者约为2 500万，在城市的十大疾病死因中居第四位，且患病率与死亡率仍呈上升趋势[2]。哮喘是一种包含气道重塑和气道感染的慢性疾病，并和气道过于紧缩和过分收缩联系在一起，我国超过2 500万人患有此病，每年死于哮喘病的人数高达60万之多[3-4]。
目前临床上主要用常规肺功能仪测定上述呼吸系统疾病，但由于这种测定需要被测试者高度配合，所以对婴幼儿、瘫痪患者、残疾人和高龄人群实施起来很困难，且检测时需有专业技师重复细致的指导才能得到准确的结果。而强迫振荡技术(forced oscillation technique, FOT)在测量中无需被测者高度配合，其基本设想是在呼吸系统中施加一定频率的振荡气流于受试者的口腔并叠加在呼吸气流之上，随气流进入气道和肺组织，测量经气道和肺组织吸收并折射后的振荡压力和振荡流量，并由此计算出一系列有用的肺阻抗参数。目前，此测量方法被证明有很好的重复性且整个过程完全无创，适合所有患者，尤其是老人、儿童和重症患者。使用单频率模式可追踪呼吸阻抗随时间的变换，而复频率模式可快速检测各个频率下的呼吸阻抗值，但国内的研究仅仅只停留在单一模式。鉴于此，本文提出一种集成单频率模式和复频率模式的强迫振荡呼吸阻抗测试仪，在此基础上增加了无线传输功能使测得的数据可以通过通用分组无线服务技术(general packet radio service，GPRS)上传到上位机，为进一步实现远程医疗奠定了基础。
1 系统方案设计
系统由下位机和上位机组成，整体设计框图如图 1所示。下位机使用功率为80 W、阻抗为4 Ω的扬声器、内径为2.2 cm的通气管路、内径为2.1 cm的高惯性管路、流量传感器和压力传感器来构建整个呼吸管路和采集系统。使用ARM STM32F4单片机为核心，利用STM32的模数转换器(digital to analog converter，DAC)产生特定频率的振荡信号，经放大电路放大后激励扬声器振动。扬声器振动箱体内的气体运动产生振荡气流，经出口和连接管进入受试者的口腔并叠加在呼吸气流之上，随气流进入气道和肺组织，由气道和肺组织吸收并折射后的压力和流量信号，经滤波电路后由STM32的模数转换器(analog to digital converter，ADC)转换成数字信号。系统在单频率模式下用移动滤波算法和互相干算法算出每个时刻对应的呼吸阻抗(resistance，Rrs)和呼吸电抗(reactance，Xrs)值并以波形的形式实时显示在触摸显示屏上，检测结束后显示整个阶段的Rrs和Xrs平均值。系统在复频率模式下用傅里叶变换将时域信号转化成频率信号，利用功率谱方法算出对应频率点的Rrs和Xrs，并在触摸显示屏上显示出频谱分析图、振荡频率为5 Hz和20 Hz下的呼吸阻抗R5和R20、振荡频率为5 Hz下的呼吸电抗X5以及共振频率Fres等。利用GPRS模块与上位机通信，建立连接后将数据上传到上位机进行进一步的处理。
图1
系统结构框图
Figure1.
System scheme
2 电路设计
系统电路主要包括STM32控制模块、振荡信号发生模块、流量和压力采集模块、外部静态随机存取存储器(static random access memory, SRAM)模块、液晶触摸显示模块、GPRS无线传输模块和电源供电模块。
2.1 STM32控制模块
本文使用ST公司基于Cortex-M4内核的32位STM32F407ZGT6单片机完成对气体流动引起的压力差电信号的采集、模数转换、数据的分析、数据的无线发送、控制对外部存储器的读写以及液晶触摸显示等功能。该单片机最高工作频率168 MHz，内带1 MB闪存，具有可变静态存储控制器(flexible static memory controller，FSMC)、定时器和3个ADC模块和2个通道DAC等外设资源。本文直接使用内置的12位DAC产生单频率和复频率信号以驱动扬声器，用内置的2个12位ADC采集气体压力和流量信号。
2.2 振荡信号发生模块
振荡信号发生模块由DAC信号发生接口、信号放大电路、重低音扬声器、圆锥箱体和气体管路组成。DAC发生的信号经放大电路后激励扬声器产生气体振荡流，推动箱体内的气流进入到管路内，再通过管路进入人体的呼吸道。
2.3 流量和压力采集模块
2.3.1 流量传感器和压力传感器的选择
流量和压力采集模块由流量传感器、压力传感器、低通滤波电路和ADC采样电路组成。流量和压力传感器将气体流量和压力转化成模拟电压信号，经二阶低通滤波电路滤除外界干扰后，由ADC采样电路采集并转化成可处理的数字信号。
使用的流量传感器是由矽翔公司生产的FS6022系列一次性呼吸传感器，可测量满量程为0~5 L/s，电压信号输出范围为0~5 V，其中0~2.5 V为正向流量，2.5~5.0 V为反向流量，电压输出与压力呈线性关系。
使用的压力传感器是由美国SMI公司生产的SM5852压差式传感器，该传感器具有多阶压力修正和温度补偿双列，确保了采集压力数据的线性和可靠性。其测量满量程是0~0.15 PSI，电压输出0.50~4.5 V，其中0.5~2.5 V为正压差，2.5~4.5 V为负压差。
2.3.2 滤波电路设计
有用的呼吸信号主要集中在1~30 Hz，因此需在流量传感器和压力传感器的输出端连接低通滤波电路以滤除高频信号干扰和50 Hz的市电干扰。本文采用的是二阶巴特沃思低通滤波电路，使用单电源电平转换芯片MAX232将+5 V的输入电压转化成10 V电压供OP07D放大器使用；由于传感器输出的电压大于ADC采样引脚可接受的电压范围(0~3.3 V)，故在其滤波电路后添加分压电路；最后得到的二阶低通滤波电路的截止频率fc=38 Hz。
2.3.3 ADC采样电路
STM32自带8通道、12位精度的ADC转换电路。ADC转换的模拟信号包括气路压力的模拟电压信号和气路流量的模拟电压信号，经二阶低通滤波电路滤波后分别由STM32外部复用端口PC0和PA0采集。采样频率由系统滴答定时器(SYSTICK)配置生成，单频率模式下振荡周期5 Hz，采样率100 Hz，复频率模式下振荡周期1~20 Hz，采样率128 Hz。
2.4 外部SRAM
STM32F407ZGT内部具有容量为192 KB的随机存取存储器(random access memory，RAM)，但由于强迫振荡需采集大小16 bit，以10 ms为间隔共32 s的流量和压力数据，在求解傅里叶变换中还需存储大量的中间变量，远超内部自带的RAM容量，故需外扩SRAM。系统选用ISSI公司生产的16位宽、1 024 K容量的静态内存芯片IS62WV102416，将该芯片连接到STM32的FSMC接口上，通过对FSMC控制器的设置可以实现在不增加外部器件的情况下对存储器的直接扩展。
2.5 液晶触摸显示模块
使用瑞佑公司生产的256/64 K色带有文字和绘图模式的双图层液晶显示驱动器RA8875，其支持16位的8080微处理机接口的传输协议。配置FSMC接口来模拟8080传输协议，选择FSMC_NE4作为片选，再加上RA8875自带触摸接口，配置完成即可通过FSMC访问RA8875相关寄存器来显示文字、图像和波形数据并可读取其触摸信号。
2.6 GPRS模块
GPRS是在现有的全球移动通信系统上发展出来的一种新的分组数据承载业务，提供端到端的连接和广域的无线IP连接。系统采用的GPRS为SIMCom公司的SIM900A模块，该模块体积小、功耗低，内嵌TCP/IP协议。STM32处理器与无线模块的物理接口为RS232，通过发送相应的“AT”指令即可完成对模块的操作。
3 算法设计
强迫振荡呼吸阻抗测量原理如图 2所示，图中人体呼吸系统总阻抗为Zrs(t)，人体与振荡源相连管路的呼吸阻抗为Z1，人体与外界相连的呼吸管路的阻抗为Z2。扬声器产生高频可控压力源Pap(tf)由振荡箱体流经呼吸管路后到达受试者内部，再提取经气道和肺组织吸收并折射的振荡压力和振荡流量信号分别为p(t, tf)和q(t, tf)。Q(t, tf)为流量传感器响应的流量，P(t, tf)为压力传感器响应的压力，根据呼吸力学模型可以得到式(1)。
图2
强迫振荡呼吸阻抗测量原理图
Figure2.
Schematic figure of forced oscillation respiratory resistance detection
| $\begin{align} & q\left( t,{{t}_{f}} \right)=\frac{1}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ & p\left( t,{{f}_{f}} \right)=\frac{Z\text{rs}}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ \end{align}$ |
由式(1)可知，虽然q(t, tf)和p(t, tf)都受振荡源和呼吸系统的影响，但其比值只受呼吸阻抗Zrs的影响，而流量传感器响应的流量和压力传感器响应的压力为人体的呼吸信号与强迫振荡信号之和，其表达式如下：
| $\begin{align} & Q\left( t,{{f}_{f}} \right)=Q\left( t \right)+q\left( t,{{t}_{f}} \right) \\ & P\left( t,{{t}_{f}} \right)=P\left( t \right)+p\left( t,{{t}_{f}} \right) \\ \end{align}$ |
为了得到呼吸阻抗Zrs，需要将呼吸信号从中分离出来。
3.1 单频率算法设计
单频率检测是在单频率振荡源下分析呼吸阻抗随时间发生的变化，信号在时域上分析处理，从而可以更直观地反映呼吸阻力在呼吸周期内发生的变化[5]。计算过程为：首先利用移动平均滤波法消除呼吸成分对呼吸阻抗计算的影响和外界的噪声干扰，再利用互相干法计算呼吸阻抗、呼吸电抗和呼吸总阻抗参数。
3.1.1 移动平均滤波
移动平均滤波技术是对有限时间内的数据信号进行求和取平均；首先取平均时间窗窗长为一个振荡周期，其对应的采样点数为N，其次对振荡周期内的采样点进行累加取平均得到ΣP(i)/N，再用该点原始采样值与该点振荡周期内的平均值进行差值得到新的数值P(i)-ΣP(i)/N。该数值即为滤除呼吸信号后的强迫振荡信号，其流程图如图 3所示。
图3
移动平均滤波法原理图
Figure3.
Schematic figure of moving average filtering method
3.1.2 互相干法
互相干法是以离散序列的离散傅里叶级数为理论基础，以其基波分量的傅里叶级数近似代替基波与各次谐波的傅里叶级数。利用三角函数形式的傅立叶展开式，因此振荡压力信号可表达为：
| $P\left( t \right)={{a}_{0}}+\sum\limits_{k=1}^{\infty }{\left[ {{a}_{k}}\cos \left( \frac{2\pi nt}{T}+{{b}_{k}}\sin \frac{2\pi nt}{T} \right) \right]}$ |
k为k阶谐波(k=1, 2, 3, …)，a0为直流分量，ak和bk为实数，与系统对外加信号的响应有关，其值由下式确定：
| $\begin{align} & {{a}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\cos \left( \frac{2\pi kt}{T} \right)\text{d}t \\ & {{b}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\sin \left( \frac{2\pi kt}{T} \right)\text{d}t \\ \end{align}$ |
运用互相干法原理，以每一时刻的基波分量(k=1)的傅里叶级数近似替代基波与各次谐波的傅里叶级数。设一个振荡周期T内有N个采样点，故在t时刻有[6]：
| $\begin{align} & {{A}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\sin \left( \frac{2\pi t}{N} \right) \\ & {{B}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\cos \left( \frac{2\pi t}{N} \right) \\ \end{align}$ |
在t时刻的振荡信号振幅和相位为：
| $\begin{align} & {{\left| p \right|}_{t}}=\sqrt{A_{t}^{2}+B_{t}^{2}} \\ & \theta \left( {{p}_{t}} \right)=\arctan \left( \frac{{{B}_{t}}}{{{A}_{t}}} \right) \\ \end{align}$ |
取移动平均滤波后的采样点t，并以该点为中心将前后半个振荡周期内的采样点数内每点乘以sin(2πt/N)并累加，再除以2/N得到At，按照同样的思路可以得到该点的Bt，At和Bt的平方并相加取平方根后得到该点的幅值｜p｜t，由arctan(Bt/At)可以计算出该点的相位值，其互相干法流程示意图如图 4所示。
图4
互相干法原理图
Figure4.
Schematic figure of cross-correlation
按照同样的方法可以计算出流量振幅｜q｜t和流量相位θ(qt)。再按照公式(7)可计算出相应呼吸总阻抗和相位信息。
| $\begin{align} & {{\left| Z \right|}_{t}}=\frac{{{\left| p \right|}_{t}}}{{{\left| q \right|}_{t}}} \\ & {{\Phi }_{t}}=\theta \left( {{p}_{t}} \right)-\theta \left( {{q}_{t}} \right) \\ \end{align}$ |
由极坐标公式可以得到：
| $\begin{align} & R\text{rs=}{{Z}_{t}}*\cos \left( {{\Phi }_{t}} \right) \\ & X\text{rs=}{{Z}_{t}}*\sin \left( {{\Phi }_{t}} \right) \\ \end{align}$ |
依照上述公式可计算出每个时刻的Rrs和Xrs。
3.2 复频率算法和设计
频谱分析法即以横坐标表示频率，纵坐标表示对应频率的函数幅值，可简便描述各频率下的频谱特性。而从时域到频域，需使用频谱分析技术即快速傅里叶变换，设一个离散的时域信号x(n)采样点数为N，离散傅里叶变换X(k)方程为[7]：
| $X\left( k \right)=\frac{1}{N}\sum\limits_{n=0}^{N-1}{x\left( n \right){{e}^{-j2\pi nt/N}}=A\left( k \right)+jB\left( k \right)}$ |
式(9)中，X(k)为频域值，k为频率值的索引k=(0, …, N 1)，A(k)为傅里叶变换的实部，B(k)为傅里叶变换的虚部。将采集到的数据进行基4快速傅里叶变换，则得到由各个频率量的实部和虚部组成的数组，并从数组中按其对应下标序列来提取出各个频率量的实部和虚部。由于振荡成分中混入一定的呼吸高频成分，为避免呼吸高频成分的影响，本文采用近似替代法，即用该频率点的往左和往右的第二位信号相加取平均来近似替代在该频率量的呼吸高频成分。用频谱中的该点减去该呼吸高频成分，便可提取出振荡波成分，其原理图如图 5所示。
图5
提取振荡成分原理图
Figure5.
Schematic figure of extracting oscillatory component
设提取到的压力和流量的频谱实部分别为Ap(fk)、Aq(fk)，虚部分别为Bp(fk)、Bq(fk)，则可得到其自功率谱和互功率谱：
| $\begin{align} & {{G}_{\text{pp}}}=\left( {{A}_{\text{p}}}+\text{j}{{B}_{\text{p}}} \right)*\left( {{A}_{\text{p}}}-\text{j}{{B}_{\text{p}}} \right)=A_{\text{p}}^{2}+B_{\text{p}}^{2} \\ & {{G}_{\text{qq}}}=\left( {{A}_{\text{q}}}+\text{j}{{B}_{\text{q}}} \right)*\left( {{A}_{\text{q}}}-\text{j}{{B}_{\text{q}}} \right)=A_{\text{q}}^{2}+B_{\text{q}}^{2} \\ & {{G}_{\text{pq}}}=\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)+\text{j}\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right) \\ \end{align}$ |
上述式中，Gpp为呼吸压力的自功率谱；Gqq为呼吸流量的自功率谱；Gqp为呼吸压力与流量的互谱；其相应频率f下的阻抗的幅值和相位角由下列式子求出：
| $\begin{align} & {{Z}_{f}}=\frac{{{G}_{\text{pp}}}}{\left| {{G}_{\text{qp}}} \right|}=\frac{A_{\text{p}}^{2}+B_{\text{p}}^{2}}{\sqrt{{{\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)}^{2}}+{{\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right)}^{2}}}} \\ & {{\theta }_{f}}=-\arctan \left( \frac{{{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}}}{{{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}}} \right) \\ \end{align}$ |
再由下式可得到相应频率下的呼吸阻抗和呼吸电抗:
| $\begin{align} & {{R}_{f}}=\left| {{Z}_{f}} \right|\cos {{\theta }_{f}} \\ & {{X}_{f}}=\left| {{Z}_{f}} \right|\sin {{\theta }_{f}} \\ \end{align}$ |
由于测量呼吸系统压力和流量信号时存在自主呼吸源等无关的噪声干扰，因此还应当引入频域相关函数r2(Coherence function)来计算其准确性
| ${{r}^{2}}=\frac{{{\left| {{G}_{\text{qp}}} \right|}^{2}}}{{{G}_{\text{pp}}}{{G}_{\text{qq}}}}$ |
r2是表征线性系统输入、输出关系的指标，r2值介于0~1之间，通常r2 > 0.95时所测得参数值才是可靠的。
4 系统性能测试
系统实物如图 6所示，将管路连接到呼吸管路接口后在显示屏上触摸相应的测试选项即可进行测试。为了评价整个系统的可靠性还需进行一些测试验证，包括机械阻抗模拟测试、模拟肺重复性测试以及初步在体测试。
图6
系统实物图
Figure6.
Photo of the system
4.1 机械阻抗测试
为了验证系统的可靠性，将塑料软管连接至呼吸管路接口上测定，根据空气动力学原理，气流流过圆管时的流体雷诺数(Reynolds number，Re)方程为[8]：
| $Re=\frac{\rho dQ}{\eta S}*{{10}^{3}}$ |
式中d：塑料管直径(cm)，Q流量(L/s)，S：管截面积(cm2)，ρ：气体密度(g/cm3)，η：气体粘滞系数(g/cm)。在测试中我们使用的流量不超过0.2 L/s，选用的管直径d=1.5 cm，取ρ空气=1.29×10-3 g/cm3，η=1.81×10-4 g/cm代入上式，得雷诺系数Re < 1997 < Re(K)(对圆管，临界雷诺数Re(K)一般取2 300)，因此在层流范围内。由公式L=ρl/r2，R=8ηL/πr4(其中：R为管阻力，L为管惯性，η为粘滞系数，r为管半径，l为管长度)，可以计算阻抗。但Z与气体成分、温度以及管子连接等机械结构有关，绝对值难以与实测值精确比较，而对与同一根塑料管，Z正比于l，故可依此做相对检验，测量的不同长度对应的Z值如表 1所示。
表1
不同长度对应的Z值
Table1.
Corresponding Z values of different lengths
用Matlab软件对呼吸阻抗测试仪测得的Z值对长度l做回归得：Z=1.093+0.049 4*l, r=0.993。可见阻抗值Z与阻抗管长度成正比，两者的线性相关性为r=0.993。
4.2 用模拟肺和志愿者测试可靠性
为了验证系统的重复性，将呼吸机和模拟肺连接至系统，设置呼吸机为CPAP模式，压力6 cm H2O，呼吸频率为12次/min，用复频率方式测试其32 s内的平均阻抗值，重复5次，得到测量值如表 2所示。用统计软件SPSS 19.0对数据进行相关分析，得到的参数用均值±标准差(x±s)形式，可见基本上满足一致性要求。
表2
用复频率方式测试模拟肺呼吸阻抗参数
Table2.
Resistance parameters of simulated lung respiration obtained using complex frequency method
为了验证系统的有效性，征集2组成年人，第一组为哮喘病患者5例，第二组为正常人5例，分别在本文研制的系统上用单频率模式下测试3次，数据结果采用均值±标准差形式，如图 7所示。
图7
5例哮喘患者和5例正常人的测试结果(x±s)
Figure7.
Test results of 5 patients with asthma and 5 healthy people (x±s)
从图中可以看出：哮喘患者的Rrs和Xrs明显高于正常人，对两组数据进行配对t检验，均得到P < 0.01，表明两组数据差异有统计学意义。
5 总结
本文以STM32单片机为主控芯片，控制呼吸流量信号的采集、处理、显示和无线传输。实验结果表明仪器不仅可以使用单频率模式追踪呼吸阻抗随时间的变化，也可以使用复频率模式快速检测各个频率上的呼吸阻抗并显示出频谱分析图。仪器提供的UI界面可为用户提供可视化操作，且测量的肺阻抗重复性高、可靠性好，具有监测肺部状况的意义。此外该仪器还具有远程传输的功能，可为进一步研究呼吸阻抗测试仪在家用市场上的发展提供参考。
引言
慢性阻塞性肺疾病(chronic obstructive pulmonary disease，COPD)和哮喘(asthma)等是常见的呼吸系统疾病。COPD是一种慢性肺部疾病，主要表现为呼出气流受限不完全可逆和呼吸困难[1]。在我国，COPD在15岁以上的人群患病率为3%，全国患者约为2 500万，在城市的十大疾病死因中居第四位，且患病率与死亡率仍呈上升趋势[2]。哮喘是一种包含气道重塑和气道感染的慢性疾病，并和气道过于紧缩和过分收缩联系在一起，我国超过2 500万人患有此病，每年死于哮喘病的人数高达60万之多[3-4]。
目前临床上主要用常规肺功能仪测定上述呼吸系统疾病，但由于这种测定需要被测试者高度配合，所以对婴幼儿、瘫痪患者、残疾人和高龄人群实施起来很困难，且检测时需有专业技师重复细致的指导才能得到准确的结果。而强迫振荡技术(forced oscillation technique, FOT)在测量中无需被测者高度配合，其基本设想是在呼吸系统中施加一定频率的振荡气流于受试者的口腔并叠加在呼吸气流之上，随气流进入气道和肺组织，测量经气道和肺组织吸收并折射后的振荡压力和振荡流量，并由此计算出一系列有用的肺阻抗参数。目前，此测量方法被证明有很好的重复性且整个过程完全无创，适合所有患者，尤其是老人、儿童和重症患者。使用单频率模式可追踪呼吸阻抗随时间的变换，而复频率模式可快速检测各个频率下的呼吸阻抗值，但国内的研究仅仅只停留在单一模式。鉴于此，本文提出一种集成单频率模式和复频率模式的强迫振荡呼吸阻抗测试仪，在此基础上增加了无线传输功能使测得的数据可以通过通用分组无线服务技术(general packet radio service，GPRS)上传到上位机，为进一步实现远程医疗奠定了基础。
1 系统方案设计
系统由下位机和上位机组成，整体设计框图如图 1所示。下位机使用功率为80 W、阻抗为4 Ω的扬声器、内径为2.2 cm的通气管路、内径为2.1 cm的高惯性管路、流量传感器和压力传感器来构建整个呼吸管路和采集系统。使用ARM STM32F4单片机为核心，利用STM32的模数转换器(digital to analog converter，DAC)产生特定频率的振荡信号，经放大电路放大后激励扬声器振动。扬声器振动箱体内的气体运动产生振荡气流，经出口和连接管进入受试者的口腔并叠加在呼吸气流之上，随气流进入气道和肺组织，由气道和肺组织吸收并折射后的压力和流量信号，经滤波电路后由STM32的模数转换器(analog to digital converter，ADC)转换成数字信号。系统在单频率模式下用移动滤波算法和互相干算法算出每个时刻对应的呼吸阻抗(resistance，Rrs)和呼吸电抗(reactance，Xrs)值并以波形的形式实时显示在触摸显示屏上，检测结束后显示整个阶段的Rrs和Xrs平均值。系统在复频率模式下用傅里叶变换将时域信号转化成频率信号，利用功率谱方法算出对应频率点的Rrs和Xrs，并在触摸显示屏上显示出频谱分析图、振荡频率为5 Hz和20 Hz下的呼吸阻抗R5和R20、振荡频率为5 Hz下的呼吸电抗X5以及共振频率Fres等。利用GPRS模块与上位机通信，建立连接后将数据上传到上位机进行进一步的处理。
图1
系统结构框图
Figure1.
System scheme
2 电路设计
系统电路主要包括STM32控制模块、振荡信号发生模块、流量和压力采集模块、外部静态随机存取存储器(static random access memory, SRAM)模块、液晶触摸显示模块、GPRS无线传输模块和电源供电模块。
2.1 STM32控制模块
本文使用ST公司基于Cortex-M4内核的32位STM32F407ZGT6单片机完成对气体流动引起的压力差电信号的采集、模数转换、数据的分析、数据的无线发送、控制对外部存储器的读写以及液晶触摸显示等功能。该单片机最高工作频率168 MHz，内带1 MB闪存，具有可变静态存储控制器(flexible static memory controller，FSMC)、定时器和3个ADC模块和2个通道DAC等外设资源。本文直接使用内置的12位DAC产生单频率和复频率信号以驱动扬声器，用内置的2个12位ADC采集气体压力和流量信号。
2.2 振荡信号发生模块
振荡信号发生模块由DAC信号发生接口、信号放大电路、重低音扬声器、圆锥箱体和气体管路组成。DAC发生的信号经放大电路后激励扬声器产生气体振荡流，推动箱体内的气流进入到管路内，再通过管路进入人体的呼吸道。
2.3 流量和压力采集模块
2.3.1 流量传感器和压力传感器的选择
流量和压力采集模块由流量传感器、压力传感器、低通滤波电路和ADC采样电路组成。流量和压力传感器将气体流量和压力转化成模拟电压信号，经二阶低通滤波电路滤除外界干扰后，由ADC采样电路采集并转化成可处理的数字信号。
使用的流量传感器是由矽翔公司生产的FS6022系列一次性呼吸传感器，可测量满量程为0~5 L/s，电压信号输出范围为0~5 V，其中0~2.5 V为正向流量，2.5~5.0 V为反向流量，电压输出与压力呈线性关系。
使用的压力传感器是由美国SMI公司生产的SM5852压差式传感器，该传感器具有多阶压力修正和温度补偿双列，确保了采集压力数据的线性和可靠性。其测量满量程是0~0.15 PSI，电压输出0.50~4.5 V，其中0.5~2.5 V为正压差，2.5~4.5 V为负压差。
2.3.2 滤波电路设计
有用的呼吸信号主要集中在1~30 Hz，因此需在流量传感器和压力传感器的输出端连接低通滤波电路以滤除高频信号干扰和50 Hz的市电干扰。本文采用的是二阶巴特沃思低通滤波电路，使用单电源电平转换芯片MAX232将+5 V的输入电压转化成10 V电压供OP07D放大器使用；由于传感器输出的电压大于ADC采样引脚可接受的电压范围(0~3.3 V)，故在其滤波电路后添加分压电路；最后得到的二阶低通滤波电路的截止频率fc=38 Hz。
2.3.3 ADC采样电路
STM32自带8通道、12位精度的ADC转换电路。ADC转换的模拟信号包括气路压力的模拟电压信号和气路流量的模拟电压信号，经二阶低通滤波电路滤波后分别由STM32外部复用端口PC0和PA0采集。采样频率由系统滴答定时器(SYSTICK)配置生成，单频率模式下振荡周期5 Hz，采样率100 Hz，复频率模式下振荡周期1~20 Hz，采样率128 Hz。
2.4 外部SRAM
STM32F407ZGT内部具有容量为192 KB的随机存取存储器(random access memory，RAM)，但由于强迫振荡需采集大小16 bit，以10 ms为间隔共32 s的流量和压力数据，在求解傅里叶变换中还需存储大量的中间变量，远超内部自带的RAM容量，故需外扩SRAM。系统选用ISSI公司生产的16位宽、1 024 K容量的静态内存芯片IS62WV102416，将该芯片连接到STM32的FSMC接口上，通过对FSMC控制器的设置可以实现在不增加外部器件的情况下对存储器的直接扩展。
2.5 液晶触摸显示模块
使用瑞佑公司生产的256/64 K色带有文字和绘图模式的双图层液晶显示驱动器RA8875，其支持16位的8080微处理机接口的传输协议。配置FSMC接口来模拟8080传输协议，选择FSMC_NE4作为片选，再加上RA8875自带触摸接口，配置完成即可通过FSMC访问RA8875相关寄存器来显示文字、图像和波形数据并可读取其触摸信号。
2.6 GPRS模块
GPRS是在现有的全球移动通信系统上发展出来的一种新的分组数据承载业务，提供端到端的连接和广域的无线IP连接。系统采用的GPRS为SIMCom公司的SIM900A模块，该模块体积小、功耗低，内嵌TCP/IP协议。STM32处理器与无线模块的物理接口为RS232，通过发送相应的“AT”指令即可完成对模块的操作。
3 算法设计
强迫振荡呼吸阻抗测量原理如图 2所示，图中人体呼吸系统总阻抗为Zrs(t)，人体与振荡源相连管路的呼吸阻抗为Z1，人体与外界相连的呼吸管路的阻抗为Z2。扬声器产生高频可控压力源Pap(tf)由振荡箱体流经呼吸管路后到达受试者内部，再提取经气道和肺组织吸收并折射的振荡压力和振荡流量信号分别为p(t, tf)和q(t, tf)。Q(t, tf)为流量传感器响应的流量，P(t, tf)为压力传感器响应的压力，根据呼吸力学模型可以得到式(1)。
图2
强迫振荡呼吸阻抗测量原理图
Figure2.
Schematic figure of forced oscillation respiratory resistance detection
| $\begin{align} & q\left( t,{{t}_{f}} \right)=\frac{1}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ & p\left( t,{{f}_{f}} \right)=\frac{Z\text{rs}}{Z\text{rs+}{{Z}_{1}}}{{P}_{\text{ap}}}\left( {{t}_{f}} \right) \\ \end{align}$ |
由式(1)可知，虽然q(t, tf)和p(t, tf)都受振荡源和呼吸系统的影响，但其比值只受呼吸阻抗Zrs的影响，而流量传感器响应的流量和压力传感器响应的压力为人体的呼吸信号与强迫振荡信号之和，其表达式如下：
| $\begin{align} & Q\left( t,{{f}_{f}} \right)=Q\left( t \right)+q\left( t,{{t}_{f}} \right) \\ & P\left( t,{{t}_{f}} \right)=P\left( t \right)+p\left( t,{{t}_{f}} \right) \\ \end{align}$ |
为了得到呼吸阻抗Zrs，需要将呼吸信号从中分离出来。
3.1 单频率算法设计
单频率检测是在单频率振荡源下分析呼吸阻抗随时间发生的变化，信号在时域上分析处理，从而可以更直观地反映呼吸阻力在呼吸周期内发生的变化[5]。计算过程为：首先利用移动平均滤波法消除呼吸成分对呼吸阻抗计算的影响和外界的噪声干扰，再利用互相干法计算呼吸阻抗、呼吸电抗和呼吸总阻抗参数。
3.1.1 移动平均滤波
移动平均滤波技术是对有限时间内的数据信号进行求和取平均；首先取平均时间窗窗长为一个振荡周期，其对应的采样点数为N，其次对振荡周期内的采样点进行累加取平均得到ΣP(i)/N，再用该点原始采样值与该点振荡周期内的平均值进行差值得到新的数值P(i)-ΣP(i)/N。该数值即为滤除呼吸信号后的强迫振荡信号，其流程图如图 3所示。
图3
移动平均滤波法原理图
Figure3.
Schematic figure of moving average filtering method
3.1.2 互相干法
互相干法是以离散序列的离散傅里叶级数为理论基础，以其基波分量的傅里叶级数近似代替基波与各次谐波的傅里叶级数。利用三角函数形式的傅立叶展开式，因此振荡压力信号可表达为：
| $P\left( t \right)={{a}_{0}}+\sum\limits_{k=1}^{\infty }{\left[ {{a}_{k}}\cos \left( \frac{2\pi nt}{T}+{{b}_{k}}\sin \frac{2\pi nt}{T} \right) \right]}$ |
k为k阶谐波(k=1, 2, 3, …)，a0为直流分量，ak和bk为实数，与系统对外加信号的响应有关，其值由下式确定：
| $\begin{align} & {{a}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\cos \left( \frac{2\pi kt}{T} \right)\text{d}t \\ & {{b}_{k}}=\frac{2}{T}\int\limits_{t}^{t+T}{p\left( t \right)}\sin \left( \frac{2\pi kt}{T} \right)\text{d}t \\ \end{align}$ |
运用互相干法原理，以每一时刻的基波分量(k=1)的傅里叶级数近似替代基波与各次谐波的傅里叶级数。设一个振荡周期T内有N个采样点，故在t时刻有[6]：
| $\begin{align} & {{A}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\sin \left( \frac{2\pi t}{N} \right) \\ & {{B}_{t}}=\frac{2}{N}\int\limits_{t-\frac{2}{N}}^{t+\frac{2}{N}}{p\left( t \right)}\cos \left( \frac{2\pi t}{N} \right) \\ \end{align}$ |
在t时刻的振荡信号振幅和相位为：
| $\begin{align} & {{\left| p \right|}_{t}}=\sqrt{A_{t}^{2}+B_{t}^{2}} \\ & \theta \left( {{p}_{t}} \right)=\arctan \left( \frac{{{B}_{t}}}{{{A}_{t}}} \right) \\ \end{align}$ |
取移动平均滤波后的采样点t，并以该点为中心将前后半个振荡周期内的采样点数内每点乘以sin(2πt/N)并累加，再除以2/N得到At，按照同样的思路可以得到该点的Bt，At和Bt的平方并相加取平方根后得到该点的幅值｜p｜t，由arctan(Bt/At)可以计算出该点的相位值，其互相干法流程示意图如图 4所示。
图4
互相干法原理图
Figure4.
Schematic figure of cross-correlation
按照同样的方法可以计算出流量振幅｜q｜t和流量相位θ(qt)。再按照公式(7)可计算出相应呼吸总阻抗和相位信息。
| $\begin{align} & {{\left| Z \right|}_{t}}=\frac{{{\left| p \right|}_{t}}}{{{\left| q \right|}_{t}}} \\ & {{\Phi }_{t}}=\theta \left( {{p}_{t}} \right)-\theta \left( {{q}_{t}} \right) \\ \end{align}$ |
由极坐标公式可以得到：
| $\begin{align} & R\text{rs=}{{Z}_{t}}*\cos \left( {{\Phi }_{t}} \right) \\ & X\text{rs=}{{Z}_{t}}*\sin \left( {{\Phi }_{t}} \right) \\ \end{align}$ |
依照上述公式可计算出每个时刻的Rrs和Xrs。
3.2 复频率算法和设计
频谱分析法即以横坐标表示频率，纵坐标表示对应频率的函数幅值，可简便描述各频率下的频谱特性。而从时域到频域，需使用频谱分析技术即快速傅里叶变换，设一个离散的时域信号x(n)采样点数为N，离散傅里叶变换X(k)方程为[7]：
| $X\left( k \right)=\frac{1}{N}\sum\limits_{n=0}^{N-1}{x\left( n \right){{e}^{-j2\pi nt/N}}=A\left( k \right)+jB\left( k \right)}$ |
式(9)中，X(k)为频域值，k为频率值的索引k=(0, …, N 1)，A(k)为傅里叶变换的实部，B(k)为傅里叶变换的虚部。将采集到的数据进行基4快速傅里叶变换，则得到由各个频率量的实部和虚部组成的数组，并从数组中按其对应下标序列来提取出各个频率量的实部和虚部。由于振荡成分中混入一定的呼吸高频成分，为避免呼吸高频成分的影响，本文采用近似替代法，即用该频率点的往左和往右的第二位信号相加取平均来近似替代在该频率量的呼吸高频成分。用频谱中的该点减去该呼吸高频成分，便可提取出振荡波成分，其原理图如图 5所示。
图5
提取振荡成分原理图
Figure5.
Schematic figure of extracting oscillatory component
设提取到的压力和流量的频谱实部分别为Ap(fk)、Aq(fk)，虚部分别为Bp(fk)、Bq(fk)，则可得到其自功率谱和互功率谱：
| $\begin{align} & {{G}_{\text{pp}}}=\left( {{A}_{\text{p}}}+\text{j}{{B}_{\text{p}}} \right)*\left( {{A}_{\text{p}}}-\text{j}{{B}_{\text{p}}} \right)=A_{\text{p}}^{2}+B_{\text{p}}^{2} \\ & {{G}_{\text{qq}}}=\left( {{A}_{\text{q}}}+\text{j}{{B}_{\text{q}}} \right)*\left( {{A}_{\text{q}}}-\text{j}{{B}_{\text{q}}} \right)=A_{\text{q}}^{2}+B_{\text{q}}^{2} \\ & {{G}_{\text{pq}}}=\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)+\text{j}\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right) \\ \end{align}$ |
上述式中，Gpp为呼吸压力的自功率谱；Gqq为呼吸流量的自功率谱；Gqp为呼吸压力与流量的互谱；其相应频率f下的阻抗的幅值和相位角由下列式子求出：
| $\begin{align} & {{Z}_{f}}=\frac{{{G}_{\text{pp}}}}{\left| {{G}_{\text{qp}}} \right|}=\frac{A_{\text{p}}^{2}+B_{\text{p}}^{2}}{\sqrt{{{\left( {{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}} \right)}^{2}}+{{\left( {{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}} \right)}^{2}}}} \\ & {{\theta }_{f}}=-\arctan \left( \frac{{{A}_{\text{p}}}{{B}_{\text{q}}}-{{B}_{\text{p}}}{{A}_{\text{q}}}}{{{A}_{\text{p}}}{{A}_{\text{q}}}+{{B}_{\text{p}}}{{B}_{\text{q}}}} \right) \\ \end{align}$ |
再由下式可得到相应频率下的呼吸阻抗和呼吸电抗:
| $\begin{align} & {{R}_{f}}=\left| {{Z}_{f}} \right|\cos {{\theta }_{f}} \\ & {{X}_{f}}=\left| {{Z}_{f}} \right|\sin {{\theta }_{f}} \\ \end{align}$ |
由于测量呼吸系统压力和流量信号时存在自主呼吸源等无关的噪声干扰，因此还应当引入频域相关函数r2(Coherence function)来计算其准确性
| ${{r}^{2}}=\frac{{{\left| {{G}_{\text{qp}}} \right|}^{2}}}{{{G}_{\text{pp}}}{{G}_{\text{qq}}}}$ |
r2是表征线性系统输入、输出关系的指标，r2值介于0~1之间，通常r2 > 0.95时所测得参数值才是可靠的。
4 系统性能测试
系统实物如图 6所示，将管路连接到呼吸管路接口后在显示屏上触摸相应的测试选项即可进行测试。为了评价整个系统的可靠性还需进行一些测试验证，包括机械阻抗模拟测试、模拟肺重复性测试以及初步在体测试。
图6
系统实物图
Figure6.
Photo of the system
4.1 机械阻抗测试
为了验证系统的可靠性，将塑料软管连接至呼吸管路接口上测定，根据空气动力学原理，气流流过圆管时的流体雷诺数(Reynolds number，Re)方程为[8]：
| $Re=\frac{\rho dQ}{\eta S}*{{10}^{3}}$ |
式中d：塑料管直径(cm)，Q流量(L/s)，S：管截面积(cm2)，ρ：气体密度(g/cm3)，η：气体粘滞系数(g/cm)。在测试中我们使用的流量不超过0.2 L/s，选用的管直径d=1.5 cm，取ρ空气=1.29×10-3 g/cm3，η=1.81×10-4 g/cm代入上式，得雷诺系数Re < 1997 < Re(K)(对圆管，临界雷诺数Re(K)一般取2 300)，因此在层流范围内。由公式L=ρl/r2，R=8ηL/πr4(其中：R为管阻力，L为管惯性，η为粘滞系数，r为管半径，l为管长度)，可以计算阻抗。但Z与气体成分、温度以及管子连接等机械结构有关，绝对值难以与实测值精确比较，而对与同一根塑料管，Z正比于l，故可依此做相对检验，测量的不同长度对应的Z值如表 1所示。
表1
不同长度对应的Z值
Table1.
Corresponding Z values of different lengths
用Matlab软件对呼吸阻抗测试仪测得的Z值对长度l做回归得：Z=1.093+0.049 4*l, r=0.993。可见阻抗值Z与阻抗管长度成正比，两者的线性相关性为r=0.993。
4.2 用模拟肺和志愿者测试可靠性
为了验证系统的重复性，将呼吸机和模拟肺连接至系统，设置呼吸机为CPAP模式，压力6 cm H2O，呼吸频率为12次/min，用复频率方式测试其32 s内的平均阻抗值，重复5次，得到测量值如表 2所示。用统计软件SPSS 19.0对数据进行相关分析，得到的参数用均值±标准差(x±s)形式，可见基本上满足一致性要求。
表2
用复频率方式测试模拟肺呼吸阻抗参数
Table2.
Resistance parameters of simulated lung respiration obtained using complex frequency method
为了验证系统的有效性，征集2组成年人，第一组为哮喘病患者5例，第二组为正常人5例，分别在本文研制的系统上用单频率模式下测试3次，数据结果采用均值±标准差形式，如图 7所示。
图7
5例哮喘患者和5例正常人的测试结果(x±s)
Figure7.
Test results of 5 patients with asthma and 5 healthy people (x±s)
从图中可以看出：哮喘患者的Rrs和Xrs明显高于正常人，对两组数据进行配对t检验，均得到P < 0.01，表明两组数据差异有统计学意义。
5 总结
本文以STM32单片机为主控芯片，控制呼吸流量信号的采集、处理、显示和无线传输。实验结果表明仪器不仅可以使用单频率模式追踪呼吸阻抗随时间的变化，也可以使用复频率模式快速检测各个频率上的呼吸阻抗并显示出频谱分析图。仪器提供的UI界面可为用户提供可视化操作，且测量的肺阻抗重复性高、可靠性好，具有监测肺部状况的意义。此外该仪器还具有远程传输的功能，可为进一步研究呼吸阻抗测试仪在家用市场上的发展提供参考。
