基于STM32及DSP库对脉搏波进行傅里叶变换的方法

2021-12-02赵星

现代计算机 2021年28期

赵星

（安徽电子科学研究所，合肥 230026）

0 引言

中医是我国传统文化中的瑰宝，脉诊在中医四诊中占有重要地位，但中医对脉象的诊断主要依靠医生自身诊断经验，受到主观因素影响较大，难免会有疏漏。随着现代科学技术的快速发展，数学、物理学、电子信息学等学科在医学领域得到大量应用，以现代科技对脉搏波进行分析研究，已经是现代医学的重要课题。而当前对脉搏波的分析中大多使用小波变换、希尔伯特-黄变换等方法进行时频分析，而非使用傅里叶变换，其目的是期望突破傅里叶变换局限性，从而更加准确地描述时间和频率的关系。然而这并不意味着傅里叶变换在脉搏波分析上就毫无用武之地，在进行如脉象分辨［1］、中年慢性病［2］等方面的研究时傅里叶变换依然有自己的价值。下文就应用STM32及ARM CMSIS DSP库对脉搏波进行傅里叶变换的方法进行分析。

1 脉搏波的提取

脉搏波是由心脏搏动沿着动脉血管向外传播而形成。人体发生病变时常常会对心脏、血管、血液等产生影响，使得脉搏波的幅度、频率、特征点等发生变化，对此做研究分析就可以了解人体的相关信息。本文分析使用的脉搏波为桡动脉脉搏波，其具体提取方法为使用一个高灵敏度的压力传感器作为脉搏传感器，将传感器的敏感点对准桡动脉跳动最强处，再以绑带固定。传感器输出的差分信号为毫伏级，需经过放大、滤波等信号调理电路处理，再通过A/D转换，变为数字信号，然后才能进行相应的数学分析。整个流程如图1所示。

图1 脉搏波提取流程

图中放大器部分本文使用仪表放大器对传感器输出信号进行放大。因人体脉搏波频率成分主要分布在0～20 Hz之间，使用低通滤波器对信号滤波。信号在滤波之后，再通过一个加法器混入一个直流分量，然后进入STM32的AD部分，从而便于其进行处理。在实际应用中依据需求可使用性能更好的专用AD芯片，在此对具体电路不再赘述。

2 STM32F4系列MCU及ARM CMSIS DSP库简介

傅里叶变换的计算量随着数据量的增大而增大，通常进行研究时使用计算机进行分析计算，但当其用于小型仪器中时主芯片计算能力可能比较有限，这时就难以保证一定的实时性。STM32F4系列是意法半导体公司基于ARM Cortex-M4内核的MCU，其内置一个32位的浮点运算单元（FPU，float point unit）［3］，支持单精度浮点运算，可以被看做协处理器。因此当其开启FPU进行复杂的浮点运算时具备明显优势，其速度可超出未开启FPU时十余倍甚至更高。Cortex-M4内核具有DSP扩展，支持多种DSP指令［4］。ARM公司的CMSIS DSP库专门为其进行了优化，在进行向量处理、矩阵处理、计算三角函数、滤波函数、统计函数、傅里叶变换等方面十分高效便捷。这也是本文使用STM32F4系列及ARM CMSIS DSP库进行傅里叶变换的原因。

3 应用STM32F405RGT6及CMSIS DSP库进行傅里叶变换

STM32F405RGT6是属于STM32F4系列的MCU芯片，其工作频率高达168 MHz，具有单精度浮点单元（FPU），支持所有ARM单精度数据处理指令和数据类型。它还实现了DSP指令用于更快速地进行数字信号处理，以及一个存储器保护单元（MPU）从而增强了应用程序的安全性。意法半导体同时推出了STM32CubeIDE作为开发工具，其基于Eclipse/CDT框架和GCC工具链，并以GDB作为调试工具，本文使用的开发工具便是STM32CubeIDE。

傅里叶变换是一种研究信号成分的重要工具。在数字信号处理中经常需要用到离散傅里叶变换（DFT，discrete fourier transform）来获取信号的频域特征。但这种方法在采样点较多时计算量过大，因此通常采用快速傅里叶变换（FFT，fast fourier transform）来实现其在工程上的应用，在CMSIS DSP库也有对应的函数可以使用。

STM32F405RGT6内部具有3个12位2.4 MSPS的A/D模块，本文以此对脉搏波的模拟信号进行模数转换。脉搏波频率成分主要分布在0～20 Hz之间，根据采样定理采样频率应当高于于信号最高频率的2倍，本文使用的采样频率是200 Hz。由傅里叶变换的原理可知，对于第n个点的频率F n有：

其中F s为采样频率，N为采样点数。可见当N越大时，获得的频率分辨率越佳。这里采集1024个数据点。本文中的编程在FFT部分使用DSP库中文件“libarm_cortexM4lf_math.a”提供的相应函数，其他部分如AD采样和定时器等都使用STM32的HAL库。这里需要注意的是DSP库中文件的选择，不仅针对不同编译器的库不同，而且对同一芯片使用不同设置时选择的文件也是不同。本文测试时开启了STM32F405RGT6的FPU，选择使用“libarm_cortexM4lf_math.a”。在AD采样时首先要对AD进行初始化，完成各种基本设置，如通道、时钟源、分辨率、数据对齐方式等。然后选择一个定时器设置为5 ms一次中断，中断触发时对信号进行AD采样。以下为程序中AD采样部分的主要代码：

HAL_ADC_Start（&hadc1）; //开始AD转换

HAL_ADC_Poll ForConversion（&hadc1,50）;//等待AD转换完成

if（HAL_ADC_Get State（&AdcHandle）==HAL_ADC_STATE_EOC_REG）//检查AD转换是否完成

{

SignValue=HAL_ADC_GetValue（&hadc1）; //获取AD转换的值

}

以上代码可用在定时器中断函数、回调函数或设置标志后轮询等多种方式中，依工程实际需要设置即可。采集到的脉搏波曲线如图2所示。

图2 脉博波曲线

取得采样数据后可以根据数据质量进行一些数字滤波等处理，这里不再赘述。FFT运算函数的输入数据和输出数据使用同一块缓冲区。在进行FFT之前，还需对输入信号数据再做一次加入虚部的处理。采集到的信号为时域信号，作为FFT函数输入变量的实部，虚部为0。然后进行FFT，再把运算结果复数求模。以下为主要代码：

for（int i=0;i＜SignLength;i++） //SignLength为采样点数1024

{

ComplexSignArray［2*i］=SignArray［i］; //采集到的信号为实部

ComplexSignArray［2*i+1］=0; //虚部设为0

SignFrequency［i］=（（i+1）-1）*200/（float）Sign-Length; //计算谱线的频率

}

arm_cfft_f32（&arm_cfft_sR_f32_len1024,ComplexSignArray,0,1）; //进行FFT计算

arm_cmplx_mag_f32（ComplexSignArray,OutputArray,SignLength）; //把运算结果复数求模

上面代码中SignArray为采集数据，计算谱线频率依据公式（1）。FFT计算的函数arm_cfft_f32中第一个参数为对计算进行配置的参数，已经包含在DSP库中的头文件“arm_const_structs.h”中，直接使用即可。第二个参数为输入数据和输出数据共用的缓冲区。第三个参数为选择正向或反向转换的标志。第四个参数为输出位反转的标志。复数求模的函数arm_cmplx_mag_f32中第一个参数为输入数据，第二个参数为输出数据，第三个参数为数据长度。将基于STM32F405RGT6的DSP库进行傅里叶变换后的频谱图如图3，由于实信号傅里叶变换的对称性，为便于特征显示图中只绘制了频谱前半部分的主要谱线。

图3 基于stm32的DSP库计算的脉搏波频谱

4 应用Python的scipy库进行傅里叶变换并对比验证

对上文中基于STM32F405RGT6的DSP库的计算结果，可以使用计算机工具来进行验证。可以使用的验证工具有多种，如MATLAB，Scilab，Octave等，本文使用Python的scipy库来进行验证。

Python是一种解释性计算机脚本语言，因其拥有众多强大的第三方库常常用于数据科学研究、web开发、网络爬虫等众多领域。scipy库便是Python的一个重要的高级科学计算库，可用于信号处理、统计分析、线性代数计算、图像处理等诸多方面，可以十分方便地进行傅里叶变换计算。其傅里叶变换的函数fft（）位于fftpack子模块中。以下为主要代码：

import numpy as np #导入numpy模块

from scipy.fftpack import fft #从scipy模块的子模块fftpack中导入fft函数

import matplotlib.pyplot as plt#导入matplotlib模块的子模块pyplot

SignLength=1024 #采样点数

SignArrayTrans=fft（SignArray）#快速傅里叶变换

SignArrayAbs=np.abs（SignArrayTrans）#对快速傅里叶变换后的复数取模

SignArrayAbs/=（SignLength/2）

SignArrayAbs［0］/=2

SignFrequency=（np.arange（1,SignLength+1）-1）*200/SignLength #计算谱线的频率

np.set_printoptions（threshold=np.inf,suppress=True）

#完整打印全部变量

print（SignArrayAbs） #输出模的值