华大半导体HC32F4A0笔记（五），使用CMSIS-DSP库进行FFT运算

一、开启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)x_q * 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处出现极大值。
幅值计算公式为：A_n * 2 / N，即512 * 2 / 1024 = 1。
即得到正弦频率为50Hz，增益为1。结果与输入信号的相符。

四、实数Q15 1024FFT的运行时间

HC32F4AO主频设在200MHz，使用Systick记录Q15 1024FFT加上转换为浮点数最后算得幅值的总时间。



运行耗时在900us左右。再分开查看

FFT耗时565us，1024个Q5定点复数转成浮点型耗时85us，算幅值耗时203us。这段时间内应该发生了几次中断，所以实际的耗时应该要比它更短。

由于实际需求中不一定需要关注全频域的幅频响应，所以算幅值的203us可以得到相当的节省。