基于STM32及DSP庫(kù)對(duì)脈搏波進(jìn)行傅里葉變換的方法

趙星

（安徽電子科學(xué)研究所，合肥 230026）

0 引言

中醫(yī)是我國(guó)傳統(tǒng)文化中的瑰寶，脈診在中醫(yī)四診中占有重要地位，但中醫(yī)對(duì)脈象的診斷主要依靠醫(yī)生自身診斷經(jīng)驗(yàn)，受到主觀因素影響較大，難免會(huì)有疏漏。隨著現(xiàn)代科學(xué)技術(shù)的快速發(fā)展，數(shù)學(xué)、物理學(xué)、電子信息學(xué)等學(xué)科在醫(yī)學(xué)領(lǐng)域得到大量應(yīng)用，以現(xiàn)代科技對(duì)脈搏波進(jìn)行分析研究，已經(jīng)是現(xiàn)代醫(yī)學(xué)的重要課題。而當(dāng)前對(duì)脈搏波的分析中大多使用小波變換、希爾伯特-黃變換等方法進(jìn)行時(shí)頻分析，而非使用傅里葉變換，其目的是期望突破傅里葉變換局限性，從而更加準(zhǔn)確地描述時(shí)間和頻率的關(guān)系。然而這并不意味著傅里葉變換在脈搏波分析上就毫無(wú)用武之地，在進(jìn)行如脈象分辨［1］、中年慢性?。?］等方面的研究時(shí)傅里葉變換依然有自己的價(jià)值。下文就應(yīng)用STM32及ARM CMSIS DSP庫(kù)對(duì)脈搏波進(jìn)行傅里葉變換的方法進(jìn)行分析。

1 脈搏波的提取

脈搏波是由心臟搏動(dòng)沿著動(dòng)脈血管向外傳播而形成。人體發(fā)生病變時(shí)常常會(huì)對(duì)心臟、血管、血液等產(chǎn)生影響，使得脈搏波的幅度、頻率、特征點(diǎn)等發(fā)生變化，對(duì)此做研究分析就可以了解人體的相關(guān)信息。本文分析使用的脈搏波為橈動(dòng)脈脈搏波，其具體提取方法為使用一個(gè)高靈敏度的壓力傳感器作為脈搏傳感器，將傳感器的敏感點(diǎn)對(duì)準(zhǔn)橈動(dòng)脈跳動(dòng)最強(qiáng)處，再以綁帶固定。傳感器輸出的差分信號(hào)為毫伏級(jí)，需經(jīng)過(guò)放大、濾波等信號(hào)調(diào)理電路處理，再通過(guò)A/D轉(zhuǎn)換，變?yōu)閿?shù)字信號(hào)，然后才能進(jìn)行相應(yīng)的數(shù)學(xué)分析。整個(gè)流程如圖1所示。

圖1 脈搏波提取流程

圖中放大器部分本文使用儀表放大器對(duì)傳感器輸出信號(hào)進(jìn)行放大。因人體脈搏波頻率成分主要分布在0～20 Hz之間，使用低通濾波器對(duì)信號(hào)濾波。信號(hào)在濾波之后，再通過(guò)一個(gè)加法器混入一個(gè)直流分量，然后進(jìn)入STM32的AD部分，從而便于其進(jìn)行處理。在實(shí)際應(yīng)用中依據(jù)需求可使用性能更好的專(zhuān)用AD芯片，在此對(duì)具體電路不再贅述。

2 STM32F4系列MCU及ARM CMSIS DSP庫(kù)簡(jiǎn)介

傅里葉變換的計(jì)算量隨著數(shù)據(jù)量的增大而增大，通常進(jìn)行研究時(shí)使用計(jì)算機(jī)進(jìn)行分析計(jì)算，但當(dāng)其用于小型儀器中時(shí)主芯片計(jì)算能力可能比較有限，這時(shí)就難以保證一定的實(shí)時(shí)性。STM32F4系列是意法半導(dǎo)體公司基于ARM Cortex-M4內(nèi)核的MCU，其內(nèi)置一個(gè)32位的浮點(diǎn)運(yùn)算單元（FPU，float point unit）［3］，支持單精度浮點(diǎn)運(yùn)算，可以被看做協(xié)處理器。因此當(dāng)其開(kāi)啟FPU進(jìn)行復(fù)雜的浮點(diǎn)運(yùn)算時(shí)具備明顯優(yōu)勢(shì)，其速度可超出未開(kāi)啟FPU時(shí)十余倍甚至更高。Cortex-M4內(nèi)核具有DSP擴(kuò)展，支持多種DSP指令［4］。ARM公司的CMSIS DSP庫(kù)專(zhuān)門(mén)為其進(jìn)行了優(yōu)化，在進(jìn)行向量處理、矩陣處理、計(jì)算三角函數(shù)、濾波函數(shù)、統(tǒng)計(jì)函數(shù)、傅里葉變換等方面十分高效便捷。這也是本文使用STM32F4系列及ARM CMSIS DSP庫(kù)進(jìn)行傅里葉變換的原因。

3 應(yīng)用STM32F405RGT6及CMSIS DSP庫(kù)進(jìn)行傅里葉變換

STM32F405RGT6是屬于STM32F4系列的MCU芯片，其工作頻率高達(dá)168 MHz，具有單精度浮點(diǎn)單元（FPU），支持所有ARM單精度數(shù)據(jù)處理指令和數(shù)據(jù)類(lèi)型。它還實(shí)現(xiàn)了DSP指令用于更快速地進(jìn)行數(shù)字信號(hào)處理，以及一個(gè)存儲(chǔ)器保護(hù)單元（MPU）從而增強(qiáng)了應(yīng)用程序的安全性。意法半導(dǎo)體同時(shí)推出了STM32CubeIDE作為開(kāi)發(fā)工具，其基于Eclipse/CDT框架和GCC工具鏈，并以GDB作為調(diào)試工具，本文使用的開(kāi)發(fā)工具便是STM32CubeIDE。

傅里葉變換是一種研究信號(hào)成分的重要工具。在數(shù)字信號(hào)處理中經(jīng)常需要用到離散傅里葉變換（DFT，discrete fourier transform）來(lái)獲取信號(hào)的頻域特征。但這種方法在采樣點(diǎn)較多時(shí)計(jì)算量過(guò)大，因此通常采用快速傅里葉變換（FFT，fast fourier transform）來(lái)實(shí)現(xiàn)其在工程上的應(yīng)用，在CMSIS DSP庫(kù)也有對(duì)應(yīng)的函數(shù)可以使用。

STM32F405RGT6內(nèi)部具有3個(gè)12位2.4 MSPS的A/D模塊，本文以此對(duì)脈搏波的模擬信號(hào)進(jìn)行模數(shù)轉(zhuǎn)換。脈搏波頻率成分主要分布在0～20 Hz之間，根據(jù)采樣定理采樣頻率應(yīng)當(dāng)高于于信號(hào)最高頻率的2倍，本文使用的采樣頻率是200 Hz。由傅里葉變換的原理可知，對(duì)于第n個(gè)點(diǎn)的頻率F n有：

其中F s為采樣頻率，N為采樣點(diǎn)數(shù)?？梢?jiàn)當(dāng)N越大時(shí)，獲得的頻率分辨率越佳。這里采集1024個(gè)數(shù)據(jù)點(diǎn)。本文中的編程在FFT部分使用DSP庫(kù)中文件“l(fā)ibarm_cortexM4lf_math.a”提供的相應(yīng)函數(shù)，其他部分如AD采樣和定時(shí)器等都使用STM32的HAL庫(kù)。這里需要注意的是DSP庫(kù)中文件的選擇，不僅針對(duì)不同編譯器的庫(kù)不同，而且對(duì)同一芯片使用不同設(shè)置時(shí)選擇的文件也是不同。本文測(cè)試時(shí)開(kāi)啟了STM32F405RGT6的FPU，選擇使用“l(fā)ibarm_cortexM4lf_math.a”。在AD采樣時(shí)首先要對(duì)AD進(jìn)行初始化，完成各種基本設(shè)置，如通道、時(shí)鐘源、分辨率、數(shù)據(jù)對(duì)齊方式等。然后選擇一個(gè)定時(shí)器設(shè)置為5 ms一次中斷，中斷觸發(fā)時(shí)對(duì)信號(hào)進(jìn)行AD采樣。以下為程序中AD采樣部分的主要代碼：

HAL_ADC_Start（&hadc1）; //開(kāi)始AD轉(zhuǎn)換

HAL_ADC_Poll ForConversion（&hadc1,50）;//等 待AD轉(zhuǎn)換完成

if（HAL_ADC_Get State（&AdcHandle）==HAL_ADC_STATE_EOC_REG）//檢查AD轉(zhuǎn)換是否完成

{

SignValue=HAL_ADC_GetValue（&hadc1）; //獲取AD轉(zhuǎn)換的值

}

以上代碼可用在定時(shí)器中斷函數(shù)、回調(diào)函數(shù)或設(shè)置標(biāo)志后輪詢(xún)等多種方式中，依工程實(shí)際需要設(shè)置即可。采集到的脈搏波曲線如圖2所示。

圖2 脈博波曲線

取得采樣數(shù)據(jù)后可以根據(jù)數(shù)據(jù)質(zhì)量進(jìn)行一些數(shù)字濾波等處理，這里不再贅述。FFT運(yùn)算函數(shù)的輸入數(shù)據(jù)和輸出數(shù)據(jù)使用同一塊緩沖區(qū)。在進(jìn)行FFT之前，還需對(duì)輸入信號(hào)數(shù)據(jù)再做一次加入虛部的處理。采集到的信號(hào)為時(shí)域信號(hào)，作為FFT函數(shù)輸入變量的實(shí)部，虛部為0。然后進(jìn)行FFT，再把運(yùn)算結(jié)果復(fù)數(shù)求模。以下為主要代碼：

for（int i=0;i＜SignLength;i++） //SignLength為采樣點(diǎn)數(shù)1024

{

ComplexSignArray［2*i］=SignArray［i］; //采 集 到的信號(hào)為實(shí)部

ComplexSignArray［2*i+1］=0; //虛部設(shè)為0

SignFrequency［i］=（（i+1）-1）*200/（float）Sign-Length; //計(jì)算譜線的頻率

}

arm_cfft_f32（&arm_cfft_sR_f32_len1024,ComplexSignArray,0,1）; //進(jìn)行FFT計(jì)算

arm_cmplx_mag_f32（ComplexSignArray,OutputArray,SignLength）; //把運(yùn)算結(jié)果復(fù)數(shù)求模

上面代碼中SignArray為采集數(shù)據(jù)，計(jì)算譜線頻率依據(jù)公式（1）。FFT計(jì)算的函數(shù)arm_cfft_f32中第一個(gè)參數(shù)為對(duì)計(jì)算進(jìn)行配置的參數(shù)，已經(jīng)包含在DSP庫(kù)中的頭文件“arm_const_structs.h”中，直接使用即可。第二個(gè)參數(shù)為輸入數(shù)據(jù)和輸出數(shù)據(jù)共用的緩沖區(qū)。第三個(gè)參數(shù)為選擇正向或反向轉(zhuǎn)換的標(biāo)志。第四個(gè)參數(shù)為輸出位反轉(zhuǎn)的標(biāo)志。復(fù)數(shù)求模的函數(shù)arm_cmplx_mag_f32中第一個(gè)參數(shù)為輸入數(shù)據(jù)，第二個(gè)參數(shù)為輸出數(shù)據(jù)，第三個(gè)參數(shù)為數(shù)據(jù)長(zhǎng)度。將基于STM32F405RGT6的DSP庫(kù)進(jìn)行傅里葉變換后的頻譜圖如圖3，由于實(shí)信號(hào)傅里葉變換的對(duì)稱(chēng)性，為便于特征顯示圖中只繪制了頻譜前半部分的主要譜線。

圖3 基于stm32的DSP庫(kù)計(jì)算的脈搏波頻譜

4 應(yīng)用Python的scipy庫(kù)進(jìn)行傅里葉變換并對(duì)比驗(yàn)證

對(duì)上文中基于STM32F405RGT6的DSP庫(kù)的計(jì)算結(jié)果，可以使用計(jì)算機(jī)工具來(lái)進(jìn)行驗(yàn)證?？梢允褂玫尿?yàn)證工具有多種，如MATLAB，Scilab，Octave等，本文使用Python的scipy庫(kù)來(lái)進(jìn)行驗(yàn)證。

Python是一種解釋性計(jì)算機(jī)腳本語(yǔ)言，因其擁有眾多強(qiáng)大的第三方庫(kù)常常用于數(shù)據(jù)科學(xué)研究、web開(kāi)發(fā)、網(wǎng)絡(luò)爬蟲(chóng)等眾多領(lǐng)域。scipy庫(kù)便是Python的一個(gè)重要的高級(jí)科學(xué)計(jì)算庫(kù)，可用于信號(hào)處理、統(tǒng)計(jì)分析、線性代數(shù)計(jì)算、圖像處理等諸多方面，可以十分方便地進(jìn)行傅里葉變換計(jì)算。其傅里葉變換的函數(shù)fft（）位于fftpack子模塊中。以下為主要代碼：

import numpy as np #導(dǎo)入numpy模塊

from scipy.fftpack import fft #從scipy模塊的子模塊fftpack中導(dǎo)入fft函數(shù)

import matplotlib.pyplot as plt#導(dǎo)入matplotlib模塊的子模塊pyplot

SignLength=1024 #采樣點(diǎn)數(shù)

SignArrayTrans=fft（SignArray）#快速傅里葉變換

SignArrayAbs=np.abs（SignArrayTrans）#對(duì)快速傅里葉變換后的復(fù)數(shù)取模

SignArrayAbs/=（SignLength/2）

SignArrayAbs［0］/=2

SignFrequency=（np.arange（1,SignLength+1）-1）*200/SignLength #計(jì)算譜線的頻率

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

#完整打印全部變量

print（SignArrayAbs） #輸出模的值

print（SignFrequency）#輸出譜線頻率

plt.plot（SignFrequency［1:（int）（SignLength/2）］,SignArrayAbs［1:（int）（SignLength/2）］,label=′脈搏波頻譜′）#繪制頻譜圖

plt.rcParams［′font.sans-serif′］=［′SimHei′］

plt.xlabel（′頻率/Hz′）

plt.ylabel（′幅值′）

plt.legend（）

plt.show（）

上面代碼中SignArray為采樣數(shù)據(jù)儲(chǔ)存的數(shù)組，譜線的頻率根據(jù)公式（1）計(jì)算。下面圖4是基于Python的scipy庫(kù)計(jì)算的脈搏波頻譜。由于實(shí)信號(hào)傅里葉變換的對(duì)稱(chēng)性，為便于特征顯示圖中只繪制了頻譜前半部分的主要譜線。

圖4 基于Python的scipy庫(kù)計(jì)算的脈搏波頻譜

將圖3與圖4對(duì)比可見(jiàn)兩圖一致。將基于STM32F405RGT6的DSP庫(kù)和Python的scipy庫(kù)分別計(jì)算的部分主要譜峰幅值及對(duì)應(yīng)頻率進(jìn)行對(duì)比，如表1所示。

表1 主要譜峰幅值及對(duì)應(yīng)頻率進(jìn)行對(duì)比

由表1可見(jiàn)兩者數(shù)據(jù)相符，其誤差已經(jīng)非常小，在應(yīng)用中可以忽略不計(jì)。

5 結(jié)語(yǔ)

本文介紹了基于STM32對(duì)脈搏波進(jìn)行采集處理及應(yīng)用DSP庫(kù)進(jìn)行傅里葉變換的方法，并將處理的結(jié)果與使用Python的scipy庫(kù)進(jìn)行傅里葉變換的結(jié)果進(jìn)行對(duì)比，驗(yàn)證了此方法的可行性，其可用于醫(yī)療儀器對(duì)脈搏波進(jìn)行傅里葉變換的分析，從而提取出有價(jià)值的人體特征參數(shù)，其處理過(guò)程對(duì)類(lèi)似的應(yīng)用也具有參考意義。