基于STM32及DSP庫對脈搏波進行傅里葉變換的方法

趙星

（安徽電子科學研究所，合肥 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） #輸出模的值

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為采樣數據儲存的數組，譜線的頻率根據公式（1）計算。下面圖4是基于Python的scipy庫計算的脈搏波頻譜。由于實信號傅里葉變換的對稱性，為便于特征顯示圖中只繪制了頻譜前半部分的主要譜線。

圖4 基于Python的scipy庫計算的脈搏波頻譜

將圖3與圖4對比可見兩圖一致。將基于STM32F405RGT6的DSP庫和Python的scipy庫分別計算的部分主要譜峰幅值及對應頻率進行對比，如表1所示。

表1 主要譜峰幅值及對應頻率進行對比

由表1可見兩者數據相符，其誤差已經非常小，在應用中可以忽略不計。

5 結語

本文介紹了基于STM32對脈搏波進行采集處理及應用DSP庫進行傅里葉變換的方法，并將處理的結果與使用Python的scipy庫進行傅里葉變換的結果進行對比，驗證了此方法的可行性，其可用于醫療儀器對脈搏波進行傅里葉變換的分析，從而提取出有價值的人體特征參數，其處理過程對類似的應用也具有參考意義。