如何使用STM32提供的DSP库函数进行FFT

　　在STM32F1系列单片机上面实现FFT
　　最近需要做一台基于stm32的示波器。如果测量信号参数，用单片机上面一些传统的方法局限性还是比较大，就开始琢磨直接上FFT。本文将以一个实例来介绍如何使用STM32提供的DSP库函数进行FFT。
　　开始是下载stm32的DSP库，提供一个下载地址：
　　下载并解压文件：
　　
　　可以看到官方提供的是256点和1024点FFT，关于FFT的点数，我以前有个误区，以为点数和采样频率有关，看了相关资料后发现两者没有本质上的关系，但是点数影响测量的分辨率。
　　将DSP库移植到Keil工程里面：
　　
　　分析一下两个FFT函数：
　　void cr4_fft_256_stm32（void *pssOUT， void *pssIN， u16 Nbin）;
　　void cr4_fft_1024_stm32（void *pssOUT， void *pssIN， u16 Nbin）;
　　// *pssOUT是FFT之后输出频域的数组，*pssIN为输入的时域采样信号数组，Nbin为FFT点数
　　在信号与系统或数字信号处理课程里面经常遇到FFT相关运算题目，根据经验，时域信号通常是实数，但是计算后得到的频域数据往往是虚数，同样按照官方库的说明，输入、输出两个数组都应该是32位的，高16位存储实部，低16位存储虚部。因此在定义数组时：
　　// long是32位
　　long InBufArray［NPT］={0}; //定义输入数组
　　long OutBufArray［NPT/2］; //定义输出数组
　　long MagBufArray［NPT/2］; //幅值
　　此库只能单纯的进行FFT运算，如果想要得到信号的频率等信息还要计算各次谐波的幅度值。需用到以下函数：
　　void GetPowerMag（）
　　{
　　signed short lX，lY;
　　float X，Y，Mag;
　　unsigned short i;
　　for（i=0; i《NPT/2; i++）
　　{
　　lX = （OutBufArray［i］ 《《 16） 》》 16;
　　lY = （OutBufArray［i］ 》》 16）;
　　//除以32768再乘65536是为了符合浮点数计算规律
　　X = NPT * （（float）lX） / 32768;
　　Y = NPT * （（float）lY） / 32768;
　　Mag = sqrt（X * X + Y * Y） / NPT;
　　if（i == 0）
　　MagBufArray［i］ = （unsigned long）（Mag * 32768）;
　　else
　　MagBufArray［i］ = （unsigned long）（Mag * 65536）; //MagBufArray［］为计算输出的幅值数组
　　}
　　}
　　有了上面的铺垫后就可以用单片机进行数据采集，然后用FFT进行处理。本实验用的单片机是普中科技stm32f103zet6的开发板（带彩屏）、函数发生器一台。
　　用stm32f103自带的12位ADC进行数据采集，FFT之后如果要获取信号频率、幅度等信息还要知道采样频率Fs，因此一般都是用定时器触发ADC采集，再用DMA进行搬运，定时器时间可以自己设置，采样频率也就知道了。具体配置如下（只贴出关键部分）：
　　//定时器1初始化
　　void TIM1_Init（void）
　　{
　　TIM_TimeBaseStructure.TIM_Period = ARR; //用了宏定义
　　TIM_TimeBaseStructure.TIM_Prescaler = PSC;
　　TIM_TimeBaseStructure.TIM_ClockDivision = 0x00;
　　TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;
　　TIM_TimeBaseInit（TIM1， &TIM_TimeBaseStructure）;
　　TIM_OCInitStructure.TIM_OCMode = TIM_OCMode_PWM1; //要配成PWM模式，用来触发ADC1
　　TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;
　　TIM_OCInitStructure.TIM_Pulse = 60;
　　TIM_OCInitStructure.TIM_OCPolarity = TIM_OCPolarity_Low;
　　TIM_OC1Init（TIM1， &TIM_OCInitStructure）;
　　TIM_CtrlPWMOutputs（TIM1， ENABLE）; //使能
　　TIM_Cmd（TIM1， DISABLE）;
　　}
　　注意定时器1在配置完成后不能马上使能，而要失能，在所有初始化完成后再使能。
　　/*********ADC1的配置********/
　　void ADC1_Init（void）
　　{
　　ADC_InitStructure.ADC_Mode = ADC_Mode_Independent;
　　ADC_InitStructure.ADC_ScanConvMode = DISABLE;
　　ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
　　ADC_InitStructure.ADC_ExternalTrigConv = ADC_ExternalTrigConv_T1_CC1; //由定时器1触发
　　ADC_InitStructure.ADC_DataAlign = ADC_DataAlign_Right;
　　ADC_InitStructure.ADC_NbrOfChannel = 1;
　　ADC_Init（ADC1， &ADC_InitStructure）;
　　ADC_RegularChannelConfig（ADC1， ADC_Channel_11， 1， ADC_SampleTime_1Cycles5）; //采样周期取最短
　　ADC_Cmd（ADC1， ENABLE）;
　　ADC_ExternalTrigConvCmd（ADC1， ENABLE）; //外部触发使能
　　ADC_DMACmd（ADC1， ENABLE）;
　　}
　　void ADC1_DMA1_Init（void）
　　{
　　NVIC_PriorityGroupConfig（NVIC_PriorityGroup_1）;
　　NVIC_InitStructure.NVIC_IRQChannel = DMA1_Channel1_IRQn; //DMA中断服务函数
　　NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 0;
　　NVIC_InitStructure.NVIC_IRQChannelSubPriority = 0;
　　NVIC_InitStructure.NVIC_IRQChannelCmd = ENABLE;
　　NVIC_Init（&NVIC_InitStructure）;
　　DMA_DeInit（DMA1_Channel1）;
　　DMA_InitStructure.DMA_PeripheralBaseAddr = ADC1_DR_Address; //ADC1对应地址
　　DMA_InitStructure.DMA_MemoryBaseAddr = （uint32_t）ADC_Value; //存储数据数组地址
　　DMA_InitStructure.DMA_DIR = DMA_DIR_PeripheralSRC;
　　DMA_InitStructure.DMA_BufferSize = ADC_LEN; //数据长度
　　DMA_InitStructure.DMA_PeripheralInc = DMA_PeripheralInc_Disable;
　　DMA_InitStructure.DMA_MemoryInc = DMA_MemoryInc_Enable;
　　DMA_InitStructure.DMA_PeripheralDataSize = DMA_PeripheralDataSize_HalfWord;
　　DMA_InitStructure.DMA_MemoryDataSize = DMA_MemoryDataSize_HalfWord;
　　DMA_InitStructure.DMA_Mode = DMA_Mode_Circular;
　　DMA_InitStructure.DMA_Priority = DMA_Priority_High;
　　DMA_InitStructure.DMA_M2M = DMA_M2M_Disable;
　　DMA_Init（DMA1_Channel1， &DMA_InitStructure）;
　　DMA_ITConfig（DMA1_Channel1， DMA_IT_TC | DMA_IT_HT， ENABLE）;
　　DMA_Cmd（DMA1_Channel1， ENABLE）;
　　}
　　此处配置了DMA中断服务函数，DMA的模式是循环读取模式，内存大小n自定义，所以单片机ADC在定时器触发之后，读取到n个数据就自动进入DMA中断服务函数。在进入服务函数之后，我们失能DMA，然后进行FFT运算，最后使能DMA。由于官方只提供了256点和1024点函数。于是我们将存储ADC数据大小的内存改为256或1024，正好可以让每次进入中断服务函数里面进行FFT的运算效果最好。
　　DMA中断服务函数：
　　void DMA1_Channel1_IRQHandler（void）
　　{
　　u16 i = 0;
　　DMA_Cmd（DMA1_Channel1， DISABLE）;
　　DMA_ClearITPendingBit（DMA_IT_HT）;
　　DMA_ClearITPendingBit（DMA1_IT_TC1）; //清除中断标志位
　　for（i=0;i《NPT;i++）
　　InBufArray［i］ = （（signed short）（ADC_Value［i］）） 《《 16; //左移16位，高16位存放实部
　　cr4_fft_1024_stm32（OutBufArray， InBufArray， NPT）;
　　GetPowerMag（）; //获取信号谐波分量的幅值
　　DMA_Cmd（DMA1_Channel1， ENABLE）;
　　}
　　由理论知识可知，幅值最大值的横坐标对应信号频率，纵坐标对应幅度。假设求得幅度值最大的数组为MagBufArray［m］=val;信号频率f0=（Fs/N）m ，信号幅值Vpp=val/（N/2）。N为FFT的点数，Fs为采样频率。如果还要求信号相位，则FFT运算后最大幅度对应的数组OutBufArray［］要转化成实部加虚部的形式：假设OutBufArray［m］=Num=a+b*i （i是虚部单位，a=Num/65536，b=Num%65536）。相位Pha=atan2（a， b），注意结果是弧度制。
　　实验验证：
　　因为短学期实习，实验室函数发生器被借走，好在之前自制了一个函数发生器（基于AD9910芯片）
　　可以通过EC11调节正弦波信号频率和幅度（幅度设置有问题，频率设置正常）。定时器的采样频率由定时器分频系数和自动装载值确定，将这两个值进行了宏定义
　　
　　假设ARR为71，PAC为9，定时器PWM频率f=72M/（（ARR+1）（PSC+1）），ADC采样频率fs=f。
　　此实验中所用采样频率为100KHz，根据奈奎斯特采样定理可知，被测信号频率最大值要小于或等于采样频率的二分之一，因此被测信号频率最大值为50Khz。用图中的41004Hz频率来测试，FFT点数为1024。运行程序，进入debug模式调试，查看数组MagBufArray最大值，可以看到：
　　
　　幅度最大值是数组的第420个，信号频率f0=（100K/1024）*420=41015Hz，精度还是很高。
　　更改采样频率为200K：
　　
　　debug：
　　
　　幅度最大值是数组的第210个，信号频率f0=（200K/1024）*210=41015Hz，误差极小。
　　更改输入信号频率为96104Hz
　　debug：
　　
　　信号频率f0=（200K/1024）*492=96093Hz，误差只有11Hz。
　　最大采样速率一直都是示波器重要的性能指标，经测试，本实验平台最大采样速率为2MHz（ARR=36-1，PSC=1-1），2M以上误差很大，被测信号最高频率为1MHz，设置输入信号
　　debug：
　　
　　f0=（2M/1024）*482=941406Hz，误差为62Hz，精度可以到0.001%。
　　分析造成误差原因：
　　FFT的分辨率，前面讲到FFT分辨率=Fs/N，要减小误差可增大N，实验用到的N=1024，Fs越大分辨率就越低，误差就越大，因此在测量低频信号时应该降低采样频率，最低要大于信号频率的两倍。
　　ADC读取数据会耗时，定时器触发AD转换时，ADC并非瞬间获取数据，而是需要有短暂的采样时间和转换时间，这样采样频率就不是理想状态下的采样频率，因此可能造成误差。
　　FFT的优缺点：
　　测量信号精度高，且最高可以测1M频率的信号，这是我上一篇文章所讲方法《用ADC和定时器测信号频率》无法比拟的，但信号频率很高，采样频率很高时，主循环将受到很大影响，因为DMA和系统总线是复用的，DMA最大能占到二分之一带宽，至少要留给系统二分之一带宽。当采样频率到2M后，由于定时器触发太频繁，导致DMA搬运ADC数据使用太多带宽，会出现主循环进不去的问题。
　　可以测量很微弱的信号，上一篇文章中讲到：信号幅度没有2V，单片机是无法检测到上升沿的，2V以下信号没办法测量频率，而FFT不一样，最低可以检测幅值为20mV的信号，包括正弦波、方波、锯齿波等等。
　　测量低频信号很耗时。