怎么实现基于STM32-AD7606的FFT交流采样？

   之前写AD7606驱动移植的时候挖下的坑，关于添加FFT计算的部分，同时本文会对AD7606的驱动程序工作方式做简要的分析。鉴于本人的水平所限，如有错误希望大家能够在评论区指正。(这次写在前面了)
   本文不会对FFT的具体原理进行讨论，有许多大佬写的非常生动形象，我会放几个我觉得不错的链接在下面。
  系统分析

  AD7606运行机制分析

   本文使用的AD7606的运行机制主要分为以下几个阶段：
   初始化阶段 开始采样阶段 结束采样阶段
   具体的各个阶段的过程由下图展现：
  
  

  

  初始化阶段主要用来设置过采样模式
  
  

  

  采样时采样频率主要由TIM4定时器来控制，定时器对时间的控制比用函数控制要精确的多，而且占用的CPU资源要少，所以一般对时间控制要求高的情况都会用定时器来控制。
  具体对TIM4频率设置的代码如下
  
void MX_TIM4_Init(uint32_t _ulFreq)
{
  TIM_ClockConfigTypeDef sClockSourceConfig = {0};
  TIM_MasterConfigTypeDef sMasterConfig = {0};


        uint16_t usPrescaler;
  uint16_t usPreiod;


  if (_ulFreq == 0)
  {
      return;                                                  //采样率为0，停止采样
  }
  else if (_ulFreq <= 100)                                     //采样频率小于100hz
  {
      usPrescaler = 42000-1;
      usPreiod = 2000 / _ulFreq;
  }
  else if(_ulFreq <= 200000)                                    //采样频率100hz - 200khz
  {
      usPrescaler = 42-1;
      usPreiod = 2000000 / _ulFreq;
  }
  else
  {
      return;
  }
       
  htim4.Instance = TIM4;
  htim4.Init.Prescaler = usPrescaler;
  htim4.Init.CounterMode = TIM_COUNTERMODE_UP;
  htim4.Init.Period = usPreiod -1;
  htim4.Init.ClockDivision = TIM_CLOCKDIVISION_DIV1;
  htim4.Init.AutoReloadPreload = TIM_AUTORELOAD_PRELOAD_ENABLE;
  if (HAL_TIM_Base_Init(&htim4) != HAL_OK)
  {
    Error_Handler();
  }
  sClockSourceConfig.ClockSource = TIM_CLOCKSOURCE_INTERNAL;
  if (HAL_TIM_ConfigClockSource(&htim4, &sClockSourceConfig) != HAL_OK)
  {
    Error_Handler();
  }
  sMasterConfig.MasterOutputTrigger = TIM_TRGO_RESET;
  sMasterConfig.MasterSlaveMode = TIM_MASTERSLAVEMODE_DISABLE;
  if (HAL_TIMEx_MasterConfigSynchronization(&htim4, &sMasterConfig) != HAL_OK)
  {
    Error_Handler();
  }


}
  前半部分是加的，后半部分使用cubemx生成的。
  这里有一个大坑需要提醒一下，单片机的IO口的速度有限，虽然AD7606只有200K的采样率，但是8通道同时采样的数据量还是比较大的，而且我用的是软件spi的模式，所以在这种应用场景下，如果以200K 8通道同时采集的模式运行的话，32会因为IO口速度的问题跳入 ERROR函数（hal库的）
  然后在定时器的中断中我们开启采集

/*
************************************************************
*函数名 ad7606_IRQSrc
*功能 定时调用函数，用于读取ADC转换数据
*形参 无
*返回值 无
************************************************************
*/
void ad7606_IRQSrc(void)
{
    uint8_t i;
    uint16_t usReadValue;
                static int j;


    /* 读取数据
    示波器监视，CS低电平持续时间 35us
    */
    AD_CS_LOW();
    for (i = 0; i < CH_NUM; i++)
    {
        usReadValue = ad7606_ReadBytes();
        if (g_tAD.usWrite < FIFO_SIZE)
        {
            g_tAD.usBuf[g_tAD.usWrite] = usReadValue;
            ++g_tAD.usWrite;
        }
    }


    AD_CS_HIGH();
//        g_tAD.usWrite = 0;
    ad7606_StartConv();
中断代码的前半部分如上，每次中断进入后会读取一个值。
  开始FFT

  数据存放

  首先我在开头定义了这些变量，用于后面计算缓存与输出，具体功能都有注释在上面了：

/*采样率*/
static uint32_t sample_freq;                     //采样率
static float32_t Freq;                            //计算出的最大频率
static float32_t Amplitude;                      //最大频率的幅值
static float32_t DC_Component;                   //直流分量


/*定义使用中的数组大小和采样点数*/
#define sample_lenth 128                         //采样点数
static float32_t InPutBuffer[sample_lenth];       //定义输入数组 InPutBuffer     该数组是外部变量
static float32_t MidBuffer[sample_lenth];         //定义中间数组
static float32_t OutPutBuffer[sample_lenth/2];    //定义输出数组 OutPutBuffer    该数组是一个静态数组
static float32_t FreqBuffer[(sample_lenth/4)-1];       //定义除了直流分量之外的频域数组


/*定义滤波中使用的中间数组*/
#define Filter_average_num 4                             //平均滤波样本数量
float32_t Mid_Filter_Freq_Buffer[Filter_average_num];    //频率平均滤波时的中间数组




/*定义fft运算中的参数*/
uint32_t fftSize = 64;                        //FFT计算点数
uint32_t ifftFlag = 0;                          //设置为fft变换   ifftFlag=1时为fft反变换
uint32_t doBitReverse = 1;                      //按位取反


/*定义计算结果存放*/
float32_t maxvalue;                             //计算出频域上最大值(包括直流分量)
uint32_t Index;                                 //最大值所在数组位置(包括直流分量)


float32_t Freq_maxvalue;                        //计算出频域上的最大值(不包括直流分量)
uint32_t Freq_Index;                            //最大值所在数组位置(不包括直流分量)


float32_t Virtual_value;                        //有效值
float32_t res;
然后我们在中断函数中加一段将采集到的数据存入数组当中。
  这里有个需要注意的地方，在dsp库中的rfft函数应用较少，cfft的函数应用较多。但是复数的表示需要用两个数组元素来表示，前一个表示实部，后面的表示虚部。
  但是对于我们采集的信号来说，我们的信号只有实部，所以我们在存取数据的时候需要将数据的虚部设置为0。
  存取代码如下：

if( j < fftSize )
                {
                       
                        InPutBuffer[2*j] = ((float)((short)g_tAD.usBuf[0])/32768/2);         //数据存放在输入数组的偶数位
                        InPutBuffer[2*j+1] = 0;                                              //奇数位置零
                        g_tAD.usWrite = 0;
                        j++;
}
我们把这一段函数放在刚才的中断函数后面。
  开始FFT计算

  这里主要用到dsp库的这几个函数
  void arm_cfft_f32(const arm_cfft_instance_f32 * S,float32_t * p1,uint8_t ifftFlag,uint8_t bitReverseFlag)
  其中 *S是用于计算的结构体，
  数组是需要计算的数组，计算的结果也会存放在这个数组中，
  ifftFlag是傅里叶变换方向的标志， 0是正向fft，1是fft逆变换
  bitReverseFlag是输出位是否翻转的标志，0为不使能，1为使能。
  用于fft运算
  void arm_cmplx_mag_f32(const float32_t * pSrc,float32_t pDst,uint32_t numSamples)*
  其中 pSrc是输入数组，
  pDst是输出数组，
  numSample是数组数量
  用于数组的求模运算
  有了以上两个函数我们就可以进行fft运算啦，想看频谱图的话只要将最后输出的数组用图片显示就可以啦。
  但是很多情况我们需要很快的知道最大的频率分量，所以我们还要用到一个函数：
  void arm_max_f32(const float32_t * pSrc,uint32_t blockSize,float32_t * pResult,uint32_t * pIndex)
  pSrc是输入数组，
  blockSize是需要计算的大小，
  pResult是数组中，
  pIndex是最大值在数组中的位置。
  这个函数可以求一个数组之中的最大值以及最大值的位置
  计算函数程序
  
arm_cfft_f32(&arm_cfft_sR_f32_len64,MidBuffer,ifftFlag,doBitReverse);      //对输入数组进行FFT变换，变换结果将存放在输入数组中
       
arm_cmplx_mag_f32(MidBuffer,OutPutBuffer,fftSize);                         //对经过FFT变换的数组进行取模运算，运算结果将存放在OutPutBuffer数组中
       
arm_max_f32(OutPutBuffer,fftSize,&maxvalue,&Index);                        //输出数组中频域最大的数值和其所在数组中的位置
频域最大的数值和其所在数组中的位置    通过这三个函数可以得到频域上幅值最大的分量的位置
  经过一系列的运算之后就可以输出啦，注释已经比较详细了，大家有不懂的可以在评论区提出
  
arm_cfft_f32(&arm_cfft_sR_f32_len64,MidBuffer,ifftFlag,doBitReverse);      //对输入数组进行FFT变换，变换结果将存放在输入数组中
       
arm_cmplx_mag_f32(MidBuffer,OutPutBuffer,fftSize);                         //对经过FFT变换的数组进行取模运算，运算结果将存放在OutPutBuffer数组中
       
arm_max_f32(OutPutBuffer,fftSize,&maxvalue,&Index);                        //输出数组中频域最大的数值和其所在数组中的位置


for(k=0;k<(fftSize/2-1);k++)
{
        FreqBuffer[k] = OutPutBuffer[k+1];                                       //取输出结果的一半，并且去除直流分量
        }
               
arm_max_f32(FreqBuffer,(fftSize/2-1),&Freq_maxvalue,&Freq_Index);          //去除直流分量后输出数组中频域最大的数值和其所在数组中的位置
               
Freq = (Freq_Index+1)*((float)sample_freq/(float)fftSize);                 //频率 = (N-1)*Fs/FFTSize        单位Hz
               
Amplitude = Freq_maxvalue/((float)fftSize/2)*10000;                        //频率幅度 = value/FFTSize/2*10   单位V
               
DC_Component = OutPutBuffer[0]/fftSize*10000;                              //直流分量 = value/FFTSize
               
               
res = ((Virtual_value-8)/43.3)/(4-((Virtual_value-8)/43.3))*2000;
               
                //printf("maxvalue = %f rn location = %d  rn",maxvalue,Index);
               
printf("Fmaxvalue = %f rn Amplitude = %f  rn  DC_Component = %f  rn  Virtual_value = %f  rn Res = %f  rn  ",Freq,Amplitude,DC_Component,Virtual_value,res);