一、开启FPU功能
点这个麻将牌四筒,展开CMSIS,把DSP勾了。
点开后
然后点这个锤子
No Auto Includes的勾不要打,让它自动include,因为CMSIS-DSP库在KEIL的安装目录中已经存在了,工程里面不需要另外添加这些库文件,自动include会帮你找到它们。
在Target标签页选上这个Use Single Precision,有些版本上面显示的是use FPU
打开hc32f4a0.h文件,搜索__FPU_PRESENT,确定其在宏定义处定义为1
至此FPU功能就开启完毕了,我们可以使用官方库函数进行多种数学和信号运算,这些官方函数在实现这些高级运算时调用了FPU的指令,这样运算速度大幅提升!
当我们使用FFT时,需要引入两个头文件,由于允许Auto Includes,KEIL MDK会自动在它的安装目录里面找到它们
#include "arm_math.h"
#include "arm_const_structs.h"
二、使用库函数q15_t arm_sin_q15(q15_t x);
生成正弦采样序列
在使用外部信号源经AD采样获取真实信号继而FFT之前,我们先模拟一个正弦波采样序列来试验测试FFT功能。
根据实际使用需求,我们模拟一个50hz的正弦信号,采样率为1600Hz,即每个周波32采样,共采样1024个点,每个点都使用Q15定点数 (-1 ~ 0.9999695)。
Fs = 1600 Hz (采样率)
N = 1024个点 (采样数)
n = 0,1,2,… N-1 (采样点序列)
t = 0,1/Fs,2/Fs,…(N-1)/Fs (采样时刻序列)
x = sin(2π x 50 x t) (采样信号序列)
#define TEST_LENGTH_SAMPLES 1024
static q15_t testInput_q15_50hz[TEST_LENGTH_SAMPLES];
uint32_t fftSize = TEST_LENGTH_SAMPLES ;
uint32_t Fs = 1600;
数学上 y = sin(x) 的定义域是是(-∞,+∞),但是DSP库中的q15_t arm_sin_q15(q15_t x);
的参数数值范围是[0,2^15),对应于弧度[0,2π)。
所以我们需要按照每周波32采样的采样率通过取余运算适配这个定义域
for(uint16_t n=0; n<fftSize; n++)
{
j = n % 32;
testInput_q15_50hz[n] = arm_sin_q15(k*j);
}
systickStart = SysTick->VAL;
其中k为:
k = 32768 * 50 / Fs;
这样我们就配置好了50Hz正弦波在1600Hz采样率下采集1024个点的信号采样序列。
三、Q15定点FFT运算
参照DSP库函数说明
/**
* @brief Instance structure for the Q15 RFFT/RIFFT function.
*/
typedef struct
{
uint32_t fftLenReal; /**< length of the real FFT. */
uint8_t ifftFlagR; /**< flag that selects forward (ifftFlagR=0) or inverse (ifftFlagR=1) transform. */
uint8_t bitReverseFlagR; /**< flag that enables (bitReverseFlagR=1) or disables (bitReverseFlagR=0) bit reversal of output. */
uint32_t twidCoefRModifier; /**< twiddle coefficient modifier that supports different size FFTs with the same twiddle factor table. */
q15_t *pTwiddleAReal; /**< points to the real twiddle factor table. */
q15_t *pTwiddleBReal; /**< points to the imag twiddle factor table. */
const arm_cfft_instance_q15 *pCfft; /**< points to the complex FFT instance. */
} arm_rfft_instance_q15;
arm_status arm_rfft_init_q15(
arm_rfft_instance_q15 * S,
uint32_t fftLenReal,
uint32_t ifftFlagR,
uint32_t bitReverseFlag);
void arm_rfft_q15(
const arm_rfft_instance_q15 * S,
q15_t * pSrc,
q15_t * pDst);
调用该函数运算FFT:
static q15_t testOutput_q15_50hz[TEST_LENGTH_SAMPLES];
arm_rfft_instance_q15 S;
uint32_t ifftFlag = 0;
uint32_t doBitReverse = 1;
/* 初始化结构体S */
(void)arm_rfft_init_q15(&S, fftSize, ifftFlag, doBitReverse);
/* 实数FFT Q15运算 */
arm_rfft_q15(&S, testInput_q15_50hz, testOutput_q15_50hz);
查阅资料得知其结果是Q5格式。
根据公式x = (float)xq * 2-Q可以换算成浮点数
static float32_t testOutput_f32_50hz[TEST_LENGTH_SAMPLES];
for(uint16_t n = 0; n < fftSize; n++)
{
testOutput_f32_50hz[n] = (float32_t)testOutput_q15_50hz[n]/32;
}
其运算结果为:
这是一个复数,偶数项为实部,其后的奇数项为虚部。这个复数为-0.03125-512j
再调用DSP库函数计算其幅值,其实由于实部很小,估算可知其幅值为512。
static float32_t testOutput[TEST_LENGTH_SAMPLES];
arm_cmplx_mag_f32(testOutput_f32_50hz, testOutput, fftSize);
求得其结果为:
由于本例是1600hz采样1024点,所以FFT后其频域响应在32 * 1600 / 1024 = 50Hz处出现极大值。
幅值计算公式为:An * 2 / N,即512 * 2 / 1024 = 1。
即得到正弦频率为50Hz,增益为1。结果与输入信号的相符。
四、实数Q15 1024FFT的运行时间
HC32F4AO主频设在200MHz,使用Systick记录Q15 1024FFT加上转换为浮点数最后算得幅值的总时间。
运行耗时在900us左右。再分开查看
FFT耗时565us,1024个Q5定点复数转成浮点型耗时85us,算幅值耗时203us。这段时间内应该发生了几次中断,所以实际的耗时应该要比它更短。
由于实际需求中不一定需要关注全频域的幅频响应,所以算幅值的203us可以得到相当的节省。