STM32H7实现全相位FFT测相（附代码）_stm32 全向傅里叶变换

一、全相位FFT的简介

1.全相位傅里叶变换（apFFT）

       传统FFT的结果可以通过一些算法实现频谱校正，如全能重心法、比值法等等。但都是基于FFT的结果，精度有限。全相位FFT算法具有初始相位不变和有效防止频谱泄露的特性。apFFT的算法是先对数据进行全相位预处理，然后进行传统的FFT运算。

2.全相位预处理过程如下：

1，构成一个N点的汉宁窗。
2，汉宁窗对自己求卷积，得到2N-1点的卷积窗。
3，求2N-1点的卷积窗的和。
4，将卷积窗的每一项除以卷积窗的和，得到2N-1点的归一化卷积窗。
5，将数据的1：2N-1项和归一化卷积窗相乘，得到加窗的2N-1项。
6，将第1项和N+1项，第2项和N+2项 ... 第N-1项和第2N-1项相加，得到经过全相预处理的N点序列。
接下来就只需要进行FFT，就可以得到apFFT的结果。

二、CubeMX配置

1.ADC同步采样：

使用ADC主从模式采样，ADC1开启DMA传输和ADC中断（ADC2不用），DMA使用normal模式，以Word字节传输，高、低16位分别存放两个adc通道数据。再以TIM3作为触发源，进行采样。

2.其他配置：

如串口重定向、DSP库配置移植可移步四臂西瓜

三、32关键代码

1.adc启动：

__IO uint8_t over_flag=0;                                            //采集完成标志

HAL_TIM_Base_Start(&htim3);  
HAL_ADCEx_Calibration_Start(&hadc1, ADC_CALIB_OFFSET_LINEARITY, ADC_SINGLE_ENDED);

2.全相位预处理

float HanNing[Length]={0};                                           //汉宁窗数组
float fft_conv[2*Length-1]={0};                                      //FFT卷积结果

void hanning_window(float *w)
{
    uint16_t n;
    for (n = 0; n < Length; n++)
    {
        w[n] = 0.5 * (1 - cos(2 * PI * n / Length));
    }
}
 
void convolve(float *w, float *y)
{
    int n, m;
    for (n = 0; n < 2 * Length - 1; n++) {
        y[n] = 0;
        for (m = 0; m < Length; m++) {
            if (n - m >= 0 && n - m < Length) {
                y[n] += w[m] * w[n - m];
            }
        }
    }
}
 
// 计算序列的和
float sum(float *x, int n)
{
    int i;
    float s = 0;
    for (i = 0; i < n; i++) {
        s += x[i];
    }
    return s;
}
 
// 归一化序列
void normalize(float *x, int n)
{
    float s = sum(x, n);
    int i;
    for (i = 0; i < n; i++) {
        x[i] /= s;
    }
}

    hanning_window(HanNing);
    convolve(HanNing,fft_conv);
    normalize(fft_conv,2*Length-1);

//相位计算函数
float Phase_atan(float32_t *inputSignal,uint32_t index)
{
    if(inputSignal[2*index+1] >= 0 && inputSignal[2*index] >= 0)
        return   0 + atan(inputSignal[2*index+1] / inputSignal[2*index]) / PI * 180;
    else if(inputSignal[2*index+1] >= 0 && inputSignal[2*index] <= 0)
        return 180 + atan(inputSignal[2*index+1] / inputSignal[2*index]) / PI * 180;
    else if(inputSignal[2*index+1] <= 0 && inputSignal[2*index] <= 0)
        return 180 + atan(inputSignal[2*index+1] / inputSignal[2*index]) / PI * 180;
    else if(inputSignal[2*index+1] <= 0 && inputSignal[2*index] >= 0)
        return 360 + atan(inputSignal[2*index+1] / inputSignal[2*index]) / PI * 180;
}

void AP_FFT(void)
{
    for (int i = 0; i < 2*Length-1; i++)
    {
        wave_2[i]=((single[i] >> 16)/65535.0f*2-1.00f)*3.57f-0.139f;
        wave_1[i]=((single[i] & 0xffff)/65535.0f*2-1.00f)*3.57f-0.139f;
    }
//    Printf(wave_1,Length);
 
    for(int i = 0; i < 2*Length-1; i++)
    {
        wave_1[i] = wave_1[i] * fft_conv[i];    //加窗
        wave_2[i] = wave_2[i] * fft_conv[i];    //加窗
    }
//   Printf(fft_conv,Length*2-1);
    for(int i = 0; i < Length-1; i++)
    {
        wave_1[i] = wave_1[i] + wave_1[Length+i];
        wave_2[i] = wave_2[i] + wave_2[Length+i];
    }
//    Printf(wave_1,Length);
    for (int i = 0; i < Length; i++)
    {
        fft_inputbuf_1[i * 2] = wave_1[i];
        fft_inputbuf_1[i * 2 + 1] = 0;
        fft_inputbuf_2[i * 2] = wave_2[i];
        fft_inputbuf_2[i * 2 + 1] = 0;
    }
    arm_cfft_radix4_f32(&scfft,fft_inputbuf_1);
    arm_cfft_radix4_f32(&scfft,fft_inputbuf_2);
 
    arm_cmplx_mag_f32(fft_inputbuf_1, fft_outputbuf_1, Length);
    arm_cmplx_mag_f32(fft_inputbuf_2, fft_outputbuf_2, Length);
 
    fft_outputbuf_1[0] = 0;
    fft_outputbuf_2[0] = 0;
 
    for(int h=0;h<Length;h++)
    {
        fft_outputbuf_1[h]/=Length/2;
        fft_outputbuf_2[h]/=Length/2;
    }
//    Printf(fft_outputbuf_1,Length);
    arm_max_f32(fft_outputbuf_1, Length/2, &Fmax_1, &Amax_pos_1);     //使用Length会让频率点在后半部分，导致频率计算错误，导致ARR->0
    arm_max_f32(fft_outputbuf_2, Length/2, &Fmax_2, &Amax_pos_2);     //使用Length会让频率点在后半部分，导致频率计算错误，导致ARR->0
    Phase_1=Phase_atan(fft_inputbuf_1,Amax_pos_1);
    Phase_1-=0.04f;                     //根据实际情况补偿，非必要不加
    Phase_2=Phase_atan(fft_inputbuf_2,Amax_pos_2);
    Phase=Phase_1-Phase_2;
//    printf("Phase1=%f Phase_2=%f Phase=%f\r\n",Phase_1,Phase_2,Phase);
    if(Phase>180) Phase = -360 + Phase;
    else if(Phase<-180) Phase = 360 + Phase;
}

主函数：

int main(void)
{
  /* USER CODE BEGIN 1 */
 
  /* USER CODE END 1 */
 
  /* MCU Configuration--------------------------------------------------------*/
 
  /* Reset of all peripherals, Initializes the Flash interface and the Systick. */
  HAL_Init();
 
  /* USER CODE BEGIN Init */
 
  /* USER CODE END Init */
 
  /* Configure the system clock */
  SystemClock_Config();
 
/* Configure the peripherals common clocks */
  PeriphCommonClock_Config();
 
  /* USER CODE BEGIN SysInit */
 
  /* USER CODE END SysInit */
 
  /* Initialize all configured peripherals */
  MX_GPIO_Init();
  MX_DMA_Init();
  MX_TIM3_Init();
  MX_USART1_UART_Init();
  MX_TIM1_Init();
  MX_DMA2D_Init();
  MX_FMC_Init();
  MX_LTDC_Init();
  MX_ADC2_Init();
  MX_ADC1_Init();
  /* USER CODE BEGIN 2 */
    RetargetInit(&huart1);
    arm_cfft_radix4_init_f32(&scfft,Length,0,1);
    HAL_TIM_Base_Start(&htim3);                                                              //1M频率
    HAL_ADCEx_Calibration_Start(&hadc1, ADC_CALIB_OFFSET_LINEARITY, ADC_SINGLE_ENDED);
    hanning_window(HanNing);
    convolve(HanNing,fft_conv);
    normalize(fft_conv,2*Length-1);
  /* USER CODE END 2 */
 
  /* Infinite loop */
  /* USER CODE BEGIN WHILE */
  while (1)
  {
      HAL_ADCEx_MultiModeStart_DMA(&hadc1,(uint32_t*)single,2*Length-1);
      while(!over_flag);
      AP_FFT();
//      printf("Phase=%f\r\n",Phase);
      over_flag=0;
 
    /* USER CODE END WHILE */
 
    /* USER CODE BEGIN 3 */
  }
  /* USER CODE END 3 */
}

四、与传统FFT对比

全相位处理：

传统FFT：

可以看出经全相位处理过后的相位非常稳定且精度大大提高。

五、结语

       本人作为新人博主，会分享自己在备战电赛时的心得体会，希望大家多多支持，评论区留企鹅号，可分享工程。

STM32F103实现双ADC同步采集电压信号

STM32ADC同步采样