开源小项目-基于STM32F1的ADC+TIM+DMA+FFT实现频率测量

2024-07-10 开源 0

一、快速傅里叶变换（FFT）

1.1 工作原理

1.2 应用

1.3 FFT官方支持库

二、使用外设简介

2.1 ADC外设

2.2 TIM外设

3.3 DMA外设

三、代码设计过程

3.1 初始化

3.2 DMA中断数据处理

3.3 其他自定义函数

1. 计算信号频率

2. 计算信号幅度

3. 计算信号相位

4. 完整代码流程

5. 注意事项

3.4 采样频率选择及配置

四、频率测量现象演示

4.1 采样频率为20Khz

4.2 采样频率为1Khz

五、FFT误差原因分析及减小误差策略

5.1 误差原因分析

5.2 FFT的优缺点补充

1. 优点

2. 缺点

该项目是基于STM32F1的频率测量，采用FFT快速傅里叶变换，配合ADC、DMA、TIM等外设实现频率测量，频率测量偏差小于千分之8。主要运用知识点：

1. ADC单通道配置与使用
2. 定时器配置及使用
3. DMA数据搬运配置与使用
4. C语言结构体、数组使用

项目源码文件下载地址：fft_measure.zip - 蓝奏云

一、快速傅里叶变换（FFT）

快速傅里叶变换（Fast Fourier Transform，简称FFT）是一种高效的算法，用于计算离散傅里叶变换（Discrete Fourier Transform，简称DFT）及其逆变换。傅里叶变换是一种将信号从时域（或空间域）转换到频域的数学方法，反之亦然。这意味着它可以将复杂的波形分解成一系列不同频率的正弦波和余弦波的组合，每个波都有其对应的振幅和相位。这一过程对于信号处理、图像处理、数据分析等领域至关重要。

1.1 工作原理

FFT通过将DFT的计算分解为较小部分的递归算法实现加速。核心思想是利用了傅里叶变换的对称性和周期性，将原问题分解为两个或更多的相似但规模较小的问题，最终通过分治策略合并这些小问题的答案以获得最终结果。特别是，它利用了蝶形运算，减少了乘法和加法的总数，这是FFT效率提升的关键。

1.2 应用

FFT在众多领域有着广泛的应用：

信号处理：分析信号的频率组成，滤波，压缩，降噪等。
图像处理：图像压缩，边缘检测，纹理分析等。
通信：频谱分析，调制解调，信号同步等。
音频处理：音乐合成，音调调整，噪声消除等。
医学：心电图、脑电图等生物信号的分析。
天文学：分析星系光谱，探测宇宙背景辐射等。

1.3 FFT官方支持库

本例中使用STM32提供的DSP库函数进行FFT运算，仅支持F1系列单片机的计算。下载链接地址：STM32_FFT官方库.zip - 蓝奏云。下载后解压文件会得到如下图1所示的源码。

其中可以看到官方提供了两个FFT计算的.s文件，FFT计算点数不同会直接影响到采集的数据精度，因此本文后面采用的是1024点的.s文件用于FFT计算。

然后打开stm32_dsp.h文件可以看到，官方提供了三个计算FFT的函数，如下：

/* Radix-4 complex FFT for STM32, in assembly  *//* 64 points*/void cr4_fft_64_stm32(void *pssOUT, void *pssIN, u16 Nbin);/* 256 points */void cr4_fft_256_stm32(void *pssOUT, void *pssIN, u16 Nbin);/* 1024 points */void cr4_fft_1024_stm32(void *pssOUT, void *pssIN, u16 Nbin);

其中每个参数的详细说明如下：

1. *pssOUT: FFT运算后输出的频域数组

2. *pssIN: 输入的时域采样信号数组

3. Nbin: 需要计算的FFT点数。

二、使用外设简介

2.1 ADC外设

模数转换，即Analog-to-Digital Converter，常称ADC，是指将连续变量的模拟信号转换为离散的数字信号的器件。其工作原理涉及对输入的模拟电压进行采样，并将其量化为一定位数的数字值。关键参数包括分辨率（通常以位数表示，如8位、10位、12位等）、转换速度（样本每秒，Sps）、精度、输入通道数量以及是否支持差分输入等。

分辨率：决定了ADC能够区分的最小电压变化量，2^n，其中n是位数。例如，一个12位的ADC可以区分2^12=4096个不同的电压等级。
转换速率（Conversion Rate/Sampling Rate）：单位时间内完成采样并转换的次数，通常以样本每秒(Samples Per Second, Sps)表示。
精度：除了分辨率外，还涉及到非线性误差、增益误差和偏置误差等因素，决定了转换结果与真实值的接近程度。
输入多路复用：许多ADC支持多个输入通道，通过多路复用器选择要转换的信号。
参考电压：决定了ADC能够测量的最大电压范围，通常是ADC的满量程电压。
电源和噪声抑制：良好的电源去耦和适当的布局对于减少噪声、提高转换精度至关重要。

其中ADC部分的主要代码初始化如下：

    /* ----------------------- ADC结构体初始化 ----------------------- */    ADC_InitTypeDef ADC_InitStructure;    ADC_InitStructure.ADC_Mode               = ADC_Mode_Independent;    ADC_InitStructure.ADC_ScanConvMode       = DISABLE;    ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;    ADC_InitStructure.ADC_ExternalTrigConv   = ADC_ExternalTrigConv_T2_CC2;  /* 使用定时器2触发ADC采样 */    ADC_InitStructure.ADC_DataAlign          = ADC_DataAlign_Right;    ADC_InitStructure.ADC_NbrOfChannel       = 1;    ADC_Init(ADC1, &ADC_InitStructure);        ADC_RegularChannelConfig(ADC1, ADC_GPIO_CH, 1, ADC_SampleTime_28Cycles5);        ADC_Cmd(ADC1,ENABLE);       /* 使能ADC */    ADC_DMACmd(ADC1, ENABLE);    ADC_ExternalTrigConvCmd(ADC1, ENABLE);  /* 外部触发使能 */

2.2 TIM外设

TIM外设，全称为Timer（定时器），是嵌入式系统中常见的外设之一，它主要用于实现定时和计数功能。TIM外设可以实现如下几种功能：

定时功能：TIM外设可以基于内部时钟或者外部输入信号进行计数，当计数值达到预设值时，可以产生中断或者DMA请求，从而触发处理器执行特定任务。这使得定时器可以用来实现精确的时间延迟、周期性任务调度等。
计数功能：连接到外部输入引脚时，TIM可以作为计数器，记录外部脉冲的数量，适用于测量频率、脉冲宽度等。
输出比较：通过比较计数器的值与预设的比较值，可以控制输出引脚的状态，实现PWM（脉宽调制）信号生成，这对于电机控制、LED亮度调节等应用非常有用。
输入捕获：捕捉外部输入信号的上升沿或下降沿，可用于测量输入信号的频率、脉冲宽度，实现编码器计数等功能。
触发和同步：高级定时器还可以与其他外设（如DAC、ADC）进行同步，用于复杂的控制应用，如在数据采集时精确控制ADC的启动时间。

其中TIM外设部分的初始化代码如下：

    /* ----------------------- TIM结构体初始化 ----------------------- */    TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure;    TIM_TimeBaseStructure.TIM_Period        = TIM_ARR - 1;          /* 自动重装载寄存器周期 */    TIM_TimeBaseStructure.TIM_Prescaler     = TIM_PSC - 1;          /* 预分频值 */    TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV1;         /* 设置时钟分割 */    TIM_TimeBaseStructure.TIM_CounterMode   = TIM_CounterMode_Up;   /* TIM向上计数模式 */    TIM_TimeBaseInit(TIM2, &TIM_TimeBaseStructure);        TIM_OCInitTypeDef TIM_OCInitStructure;    TIM_OCInitStructure.TIM_OCMode      = TIM_OCMode_PWM1;          /* 选择定时器模式:TIM脉冲宽度调制模式1 */    TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;   /* 比较输出使能 */    TIM_OCInitStructure.TIM_Pulse       = TIM_ARR/2;    TIM_OCInitStructure.TIM_OCPolarity  = TIM_OCPolarity_Low;       /* 输出极性:TIM输出比较极性低 */    TIM_OC2Init(TIM2, & TIM_OCInitStructure);                       /* 初始化外设TIM2_CH2 */        TIM_Cmd(TIM2, ENABLE);

3.3 DMA外设

DMA（Direct Memory Access，直接存储器访问）外设是一种专为高效数据传输设计的硬件控制器，它允许外部设备（如传感器、ADC、DAC、网络接口等）和系统内存（RAM）之间直接进行高速数据交换，而无需CPU的持续介入。这一机制显著减轻了CPU的负担，使得CPU可以在数据传输的同时执行其他任务，从而提高了系统的整体性能和响应速度。DMA的主要特点包括：

高效传输：DMA能够在不占用CPU资源的情况下，快速传输大量数据，特别适合于大数据量的连续传输，如音频、视频流处理，大容量存储读写等。
自动操作：DMA传输由硬件自动完成，一旦配置好DMA控制器（包括源地址、目标地址、传输量和控制参数），数据传输就会在指定的触发条件下自动开始并完成。
减少CPU干预：通过DMA，CPU不必参与每次数据的读写操作，仅在传输开始和结束时（或发生错误时）接收中断通知，这样就释放了CPU来执行其他计算密集型任务。
多通道支持：现代DMA控制器通常支持多个独立通道，每个通道可以被分配给不同的外设，同时进行数据传输，增加了系统的并行处理能力。
优先级管理：DMA控制器能够管理不同通道的优先级，确保关键数据的及时传输，特别是在多个外设同时请求DMA服务时。
灵活配置：DMA传输可以配置为单次传输或循环传输，对于需要周期性数据更新的应用（如采样数据的连续采集）非常有用。
对齐要求：为了优化传输效率，源和目标地址通常需要按照数据传输宽度（如字节、半字、全字）对齐。

其中DMA部分的主要代码初始化如下：

    /* ----------------------- DMA结构体初始化 ----------------------- */    DMA_InitTypeDef DMA_InitStructure;    DMA_InitStructure.DMA_PeripheralBaseAddr = (uint32_t)&(ADC1->DR);   /* ADC1对应地址 */    DMA_InitStructure.DMA_MemoryBaseAddr     = (uint32_t)&adcValue[0];  /* 存储数据数组地址 */    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模式 循环传输 */    DMA_InitStructure.DMA_Priority           = DMA_Priority_High;       /* 优先级 高 */    DMA_InitStructure.DMA_M2M                = DMA_M2M_Disable;         /* 禁止内存到内存的传输 */    DMA_Init(ADC_DMA_CH, &DMA_InitStructure);        DMA_ITConfig(ADC_DMA_CH, DMA_IT_TC, ENABLE); /* 使能DMA传输完成中断 */

三、代码设计过程

3.1 初始化

代码的初始化分为两部分，一部分是初始化硬件引脚信息；另一部分是初始化内置结构体信息（用于保存待测频率的基本信息）。

①硬件引脚初始化：采用STM32F1的ADC引脚（ADC1 PA1）、DMA搬运（DMA1 通道1）、TIM定时器（TIM2 通道2），整体初始化代码设计如下：（注意打开C99 Mode）

/* -------------------------- 宏定义区 -------------------------- *//* 使用的ADC引脚 PA1 */#define     ADC_GPIO_PORT   GPIOA#define     ADC_GPIO_PIN    GPIO_Pin_3#define     ADC_GPIO_CH     ADC_Channel_3/* 这里不需改动 */#define     ADC_TIM_CH      TIM_OC2Init#define     ADC_DMA_CH      DMA1_Channel1/* 定时器采集频率 */#define     TIM_ARR         (720)#define     TIM_PSC         (5)/* FFT采集点数 */#define     FFT_POINTS      (1024)#define     ADC_LEN         (FFT_POINTS)/** * @brief  adc外设引脚初始化 * @param  [in] *   @arg  None * @retval None */static void fft_adc_pin_init(void){    /* ----------------------- RCC时钟初始化 ----------------------- */    RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA, ENABLE);    RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE);    RCC_APB2PeriphClockCmd(RCC_APB2Periph_AFIO, ENABLE);    RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2, ENABLE);    RCC_AHBPeriphClockCmd(RCC_AHBPeriph_DMA1, ENABLE);    RCC_ADCCLKConfig(RCC_PCLK2_Div6);   /* 设置ADC分频因子 最大时钟不能超过14MHz */        /* ----------------------- GPIO结构体初始化 ----------------------- */    GPIO_InitTypeDef GPIO_InitStructure;    GPIO_InitStructure.GPIO_Pin   = ADC_GPIO_PIN;    GPIO_InitStructure.GPIO_Mode  = GPIO_Mode_AIN;    GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;    GPIO_Init(ADC_GPIO_PORT, &GPIO_InitStructure);        /* ----------------------- ADC结构体初始化 ----------------------- */    ADC_InitTypeDef ADC_InitStructure;    ADC_InitStructure.ADC_Mode               = ADC_Mode_Independent;    ADC_InitStructure.ADC_ScanConvMode       = DISABLE;    ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;    ADC_InitStructure.ADC_ExternalTrigConv   = ADC_ExternalTrigConv_T2_CC2;  /* 使用定时器2触发ADC采样 */    ADC_InitStructure.ADC_DataAlign          = ADC_DataAlign_Right;    ADC_InitStructure.ADC_NbrOfChannel       = 1;    ADC_Init(ADC1, &ADC_InitStructure);        ADC_RegularChannelConfig(ADC1, ADC_GPIO_CH, 1, ADC_SampleTime_28Cycles5);        ADC_Cmd(ADC1,ENABLE);       /* 使能ADC */    ADC_DMACmd(ADC1, ENABLE);    ADC_ExternalTrigConvCmd(ADC1, ENABLE);  /* 外部触发使能 */        /* ----------------------- TIM结构体初始化 ----------------------- */    TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure;    TIM_TimeBaseStructure.TIM_Period        = TIM_ARR - 1;          /* 自动重装载寄存器周期 */    TIM_TimeBaseStructure.TIM_Prescaler     = TIM_PSC - 1;          /* 预分频值 */    TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV1;         /* 设置时钟分割 */    TIM_TimeBaseStructure.TIM_CounterMode   = TIM_CounterMode_Up;   /* TIM向上计数模式 */    TIM_TimeBaseInit(TIM2, &TIM_TimeBaseStructure);        TIM_OCInitTypeDef TIM_OCInitStructure;    TIM_OCInitStructure.TIM_OCMode      = TIM_OCMode_PWM1;          /* 选择定时器模式:TIM脉冲宽度调制模式1 */    TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;   /* 比较输出使能 */    TIM_OCInitStructure.TIM_Pulse       = TIM_ARR/2;    TIM_OCInitStructure.TIM_OCPolarity  = TIM_OCPolarity_Low;       /* 输出极性:TIM输出比较极性低 */    ADC_TIM_CH(TIM2, & TIM_OCInitStructure);                        /* 初始化外设TIMx_CHx */        TIM_Cmd(TIM2, ENABLE);        /* ----------------------- DMA结构体初始化 ----------------------- */    DMA_InitTypeDef DMA_InitStructure;    DMA_InitStructure.DMA_PeripheralBaseAddr = (uint32_t)&(ADC1->DR);   /* ADC1对应地址 */    DMA_InitStructure.DMA_MemoryBaseAddr     = (uint32_t)&adcValue[0];  /* 存储数据数组地址 */    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模式 循环传输 */    DMA_InitStructure.DMA_Priority           = DMA_Priority_High;       /* 优先级 高 */    DMA_InitStructure.DMA_M2M                = DMA_M2M_Disable;         /* 禁止内存到内存的传输 */    DMA_Init(ADC_DMA_CH, &DMA_InitStructure);        DMA_ITConfig(ADC_DMA_CH, DMA_IT_TC, ENABLE); /* 使能DMA传输完成中断 */        /* ----------------------- NVIC结构体初始化 ----------------------- */    NVIC_InitTypeDef NVIC_InitStructure;    NVIC_InitStructure.NVIC_IRQChannel                   = DMA1_Channel1_IRQn;  /* DMA中断服务函数 */    NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 0;    NVIC_InitStructure.NVIC_IRQChannelSubPriority        = 3;    NVIC_InitStructure.NVIC_IRQChannelCmd                = ENABLE;    NVIC_Init(&NVIC_InitStructure);        DMA_Cmd(ADC_DMA_CH, ENABLE);        ADC_ResetCalibration(ADC1);                 /* 使能复位校准 */    while(ADC_GetResetCalibrationStatus(ADC1)); /* 等待复位校准结束 */    ADC_StartCalibration(ADC1);                 /* 开启AD校准 */    while(ADC_GetCalibrationStatus(ADC1));      /* 等待校准结束 */}

②内置结构体信息初始化：该结构体的成员变量如下：

typedef struct {    unsigned char flag;         /* DMA数据搬运完成标志位 */    unsigned short index;       /* 频率所在的fft数组索引 */    unsigned short amplitude;   /* 幅值 mV */    float sample_freq;          /* 采样频率 根据奈奎斯特采样定理 待测频率需小于等于采样频率的1/2 */    float freq;                 /* 待测频率 hz */    float phase;                /* 相位 ° */} STRUCT_WAVE_TYPEDEF;

结构体实例化代码如下：

static STRUCT_WAVE_TYPEDEF wave;                /* 定义内置结构体 保存波形信息 *//** * @brief  wave结构体初始化 * @param  [in] *   @arg  None * @retval None */static void fft_struct_init(void){    wave.flag = 0;    wave.index = FFT_POINTS / 2;    wave.sample_freq = SystemCoreClock/(TIM_ARR*TIM_PSC);        wave.freq = 500;    wave.amplitude = 3300;    wave.phase = 0;}

3.2 DMA中断数据处理

我们在通过前面的一系列配置后，完成了DMA的数据搬运配置工作，然后我们需要在DMA的中断服务函数中去拿数据，代码配置如下：

/** * @brief  DMA1通道1中断服务函数 * @param  [in] *   @arg  None * @retval None */void DMA1_Channel1_IRQHandler(void){    if(DMA_GetITStatus(DMA1_IT_TC1) == SET) {        DMA_ClearITPendingBit(DMA1_IT_TC1); /* 清除中断标志位 */        wave.flag = 1;    }}

这里主要是在中断服务函数触发中断标志后，代表ADC数据获取成功，将标志位置为1，等待后续运算。

3.3 其他自定义函数

我们通过DMA中断服务函数触发是否成功来判定ADC数据是否获取成功，当获取成功后，我们需要使用ADC的数据进行FFT的运算，该部分代码如下：

        /* 2. 将得到的ADC数据处理后放入in_buff数组中 */        for(i=0; i<FFT_POINTS; i++) {            in_buff[i] = ((signed short)(adcValue[i])) << 16;  /* 左移16位 高16位存放实部 */        }        /* 3. 调用FFT库计算 */        cr4_fft_1024_stm32(out_buff, in_buff, FFT_POINTS);

而完成FFT运算之后，如果我们需要得到信号的频率、幅值、相位等信息还要计算各次谐波的幅度值，拿到幅度值之后，就可以通过幅度值数组计算波形的信息了。该函数定义如下：

/** * @brief  获取信号谐波分量的幅值 * @param  [in] *   @arg  None * @retval None */static void fft_get_harmonic_data(void){    signed short lX, lY;    float X, Y, harmonic;    unsigned short i;        const float scale_factor = FFT_POINTS / 32768.0f;                   /* 预先计算缩放因子 避免循环内重复计算 */    const float post_scale =  (FFT_POINTS == 1) ? 32768.0f : 65536.0f;  /* 根据是否为直流分量决定后续缩放 */        for (i = 0; i < FFT_POINTS / 2; ++i) {        lX = out_buff[i] & 0xFFFF; /* 直接提取低16位 避免移位操作 */        lY = out_buff[i] >> 16;    /* 提取高16位 */        /* 使用预先计算的缩放因子进行计算 */        X = lX * scale_factor;        Y = lY * scale_factor;        harmonic = sqrtf(X * X + Y * Y); /* 使用sqrtf 针对浮点数优化的sqrt版本 */        harmonic_buff[i] = (unsigned long)(harmonic * post_scale); /* 根据索引决定最终缩放并存储 */    }    /* 特殊处理直流分量 */    if (i == 0 && FFT_POINTS > 1) {        harmonic_buff[0] /= 2;    }}

1. 计算信号频率

我们可以通过上面函数（fft_get_harmonic_data）可以获取到波形信号的所有谐波分量的幅值，当找到该幅值数组harmonic_buff中的最大幅值即为对应的待测波形频率（Fs），可通过如下公式计算得出：

Fs = (F0 / N) * max

其中F0是采样频率，N是FFT的点数，max是数组中最大幅值对应的索引。这个频率Fs即为信号的主要频率成分，即我们需要的待测频率。因此我们可以通过函数findMax_index找到数组中的最大值索引值，函数代码如下：

/** * @brief  查找数组中的最大值所在的数组索引 * @param  [in] *   @arg  *arr - 待查找数组 *   @arg  len - 查找数组长度 * @retval None */static int findMax_index(int *arr, int len){    int i;    int max_val = arr[0];    int max_index = 0;        if (len <= 0) {        return -1; /* 错误值 表示数组为空或非法大小 */    }    for (i = 1; i < len; i++) {   /* 遍历数组 */        if (arr[i] > max_val) { /* 如果找到更大的值 */            max_val = arr[i];            max_index = i;      /* 更新最大值及其索引 */        }    }    return max_index;}

2. 计算信号幅度

当我们得到信号的频率后，信号的峰值幅度可以通过harmonic_buff[max]计算得到，但通常需要根据FFT的具体实现进行适当的缩放。最终的峰峰值Vpp（即信号的最大幅度变化范围）可以通过公式获得，公式如下：

Vpp = harmonic_buff[max] / (N / 2)

这里的缩放因子N/2是因为FFT的结果包含了正负频率部分，且对于实数信号，能量分布在正负频率上是相同的，所以原始信号的能量需要平均分配到这两个部分。

3. 计算信号相位

然后我们为了获取信号的相位信息，需要利用FFT结果的复数形式。假设out_buff[max]可以表示为一个复数，拆分为实部a和虚部b，即：

out_buff[max] = a + bi

其中实部（a）和虚部（b）的计算公式如下：

实部：a = out_buff[max] >> 16

虚部：b = out_buff[max] & 0xFFFF

注意：这里假设out_buff[max]是经过适当缩放以适应整数表示的复数。因此信号的相位可以通过atan2(a, b)函数计算得到，该函数返回的是弧度值。atan2相比普通的atan函数，能正确处理所有象限的角度，避免了正负符号和0的问题。

4. 完整代码流程

/** * @brief  fft采集并计算波形信息 * @param  [in] *   @arg  None * @retval None * @note   该函数需定时执行 注意定时执行周期不宜太短 */void fft_measure_handle(void){    unsigned short i = 0;    signed short real_part, imaginary_part; /* 定义实部和虚部分量 */        if(wave.flag == 1) {        wave.flag = 0;        /* 1. 失能DMA1通道1 以防在中断过程中发生DMA传输 */        DMA_Cmd(ADC_DMA_CH, DISABLE);                /* 2. 将得到的ADC数据处理后放入in_buff数组中 */        for(i=0; i<FFT_POINTS; i++) {            in_buff[i] = ((signed short)(adcValue[i])) << 16;  /* 左移16位 高16位存放实部 */        }        /* 3. 调用FFT库计算 */        cr4_fft_1024_stm32(out_buff, in_buff, FFT_POINTS);                /* 4. 获取信号谐波分量的幅值 */        fft_get_harmonic_data();                /* 5. 获取fft频率数组最大值索引 */        wave.index = findMax_index(&harmonic_buff[1], FFT_POINTS/2-1) + 1;/* 这里有个坑 第0个数据是异常点 所以从第一个数据开始取 所以索引值index需要加1 */                /* 6. 计算波形的频率 幅值 相位等信息 */        wave.freq = wave.sample_freq / (float)FFT_POINTS * (float)(wave.index); /* 频率计算公式: Fs = F0/points*index */        wave.amplitude = harmonic_buff[wave.index] / (FFT_POINTS/2) * 70 / 100; /* 幅值计算公式: Vpp = buff[index]/(points/2) */                real_part = out_buff[wave.index] >> 16;            /* 计算实部 */        imaginary_part = out_buff[wave.index] & 0xFFFF;    /* 计算虚部 */        wave.phase = atan2(real_part, imaginary_part);     /* 相位计算公式: phase = atan2(real, imaginary) */                /* 7. 重新使能DMA1通道1 开始工作 */        DMA_Cmd(ADC_DMA_CH, ENABLE);    }}

5. 注意事项

确保在进行相位计算前，正确地将out_buff[max]转换为实部和虚部，且考虑是否需要进一步的缩放以匹配实际的物理量。
计算得到的相位可能需要根据具体应用场景调整，比如从弧度转换为度，或根据信号起始相位进行校正。
在实际应用中，还需考虑FFT窗口效果、频谱泄露等对结果的影响，可能需要采取相应的预处理或后处理措施来提高分析准确性。

3.4 采样频率选择及配置

由于在初始化中是使用的定时器+DMA采集ADC引脚输入的波形频率，因此指定定时器的频率大小就是指定采样频率的大小，而根据奈奎斯特采样定理可知，采样频率（定时器的频率）需大于等于待测频率的2倍才行，通常选择3~5倍左右。因此我们可以通过指定定时器的频率从而指定待测ADC输入波形的可被采集的频率范围。定时器频率设置如下：

/** * @brief  设置定时器采样频率(标准频率) * @param  [in] *   @arg  定时器采样频率 单位Hz * @retval None */void fft_set_sample_freq(float freq){    unsigned short psc;    wave.sample_freq = freq;    psc = SystemCoreClock/TIM_ARR/freq; /* 获取预分频值 */        if(psc < 2) { /* 预分频值小于2 采用重装载值的方式更新定时器频率 */        psc = 2;        TIM_ARRPreloadConfig(TIM2, DISABLE);        TIM2->ARR = SystemCoreClock/(psc*freq);        TIM2->CCR2 = SystemCoreClock/(psc*freq)/2;/* 保证定时器2通道2的占空比为50% */        TIM_ARRPreloadConfig(TIM2, ENABLE);    }    else { /* 预分频值大于2 采用预分频值的方式更新定时器频率 */        TIM_PrescalerConfig(TIM2, psc, TIM_PSCReloadMode_Immediate);    }}

这里需要注意：不能使用太高的定时器采样频率，因为太高的定时器采样频率会占用CPU的资源，导致定时器

四、频率测量现象演示

4.1 采样频率为20Khz

源码中定时器初始频率设置为20KHz，因此可测得的频率范围为10Khz以下，当在信号发生器一端设置波形为正弦波、频率为9500Hz、占空比为50%、幅度值为3.3V时，设置的图像见下图2：

此时使用Keil 自带的debug窗口进行观测，其测得的波形信息如下图3所示：

由上图可以得知，波形的频率与信号发生器设置的频率大体一致，由些许误差，误差率为：(9500-9492.1875) / 9500 = 0.000822，可以看到测量精度还是挺高的。

4.2 采样频率为1Khz

我们可以通过调用函数fft_set_sample_freq即可任意设置定时器的频率，即采样频率，当然也可以在初始化部分去配置定时器的重装载值和预分频值。见下图4。

接着需要在信号发生器一端设置波形：更改为方波、频率为300Hz、占空比为50%、幅度值为3.3V时，设置的图像见下图5：

接着重新进入Keil的debug进行调试，可以看到如下图6所示波形信息：

通过上图可以看到，方波的幅值和相位都比较接近信号发生器的波形数据，待测频率与信号发生器的差别也不算很大，偏差率为：(335-337.890625)/335 = -0.00863

五、FFT误差原因分析及减小误差策略

5.1 误差原因分析

FFT分辨率与采样频率：FFT的分辨率确实由采样频率F0与FFT点数N决定，即分辨率=F0/N。提高分辨率意味着能更精细地分辨相邻频率成分，在保持N不变时，增加F0会提高分辨率。而为了减小低频信号测量的误差，应当降低采样频率，但必须遵循奈奎斯特采样定理，即采样频率至少为信号最高频率的两倍。
ADC读取延迟：ADC转换存在采样时间和转换时间，这导致实际采样间隔与理论设定值不符，从而引入误差。可以通过优化ADC配置，如使用更快的采样率、减少采样保持时间、采用更高的转换精度等方法来减少这种误差。
量化误差：ADC将连续的模拟信号量化为有限位数的数字信号，此过程必然会产生量化误差，尤其是在信号幅度接近ADC分辨率极限时。选择更高位数的ADC可以减小量化误差。
频谱泄漏：非整周期信号在进行FFT时会产生频谱泄漏，导致能量分散到邻近的频率上，影响频率成分的准确识别。采用窗函数可以一定程度上减少频谱泄漏。
DC偏移和噪声：信号中的直流偏移和随机噪声也会干扰FFT分析，影响幅度和频率的准确度。预处理如高通滤波去除DC偏移、低通滤波减少噪声等可以改善分析质量。

5.2 FFT的优缺点补充

1. 优点

高精度与宽频率范围：FFT能够提供高精度的频谱分析，理论上可测量高达几百万赫兹的信号频率，远超过传统的时域测量方法，适用于复杂信号的分析。
微弱信号检测：FFT具有很高的灵敏度，即使信号幅度很小，也能通过累计多次采样来检测到，适用于微弱信号的分析，如电磁辐射、声纳信号等。
多信号同时分析：在单一FFT结果中，可以同时观察到信号中的多种频率成分，便于分析信号的组成结构。

2. 缺点

计算复杂度：虽然FFT相对于DFT大大降低了计算复杂度，但对于大数据量的处理仍然需要较大的计算资源。
低频信号处理时间长：对于低频信号，由于其周期长，需要较长的采样时间才能获得足够的周期数，导致分析过程耗时较长。
系统资源占用：在高采样率下，特别是配合DMA使用时，FFT处理和数据搬运会占用较多的系统资源，影响实时性，可能导致主循环响应变慢。
频谱泄漏与混叠：如前所述，FFT分析易受频谱泄漏和频率混叠的影响，需要合理选择窗函数和采样率来减小这些影响。

---------------------------------------------------------------------------------------------------------------------------------

分割线(以下是续写内容)

---------------------------------------------------------------------------------------------------------------------------------

前面实现了基于FFT的频率采集，后续学习了一些FFT的相关理论后，补充如下概念：

总谐波失真（Total Harmonic Distortion，简称THD）是衡量信号失真程度的重要指标，特别是在音频、电力电子设备和通信系统等领域。THD表示一个信号中谐波分量（除了基波之外的频率分量）相对于基波的比例。具体来说，THD是所有谐波分量的平方和与基波分量平方的比值的平方根，通常用百分比表示。其公式表示为：

其中：V1是基波幅度，V2,V3……Vn是谐波幅度。

void fft_measure_handle(void){    unsigned short i = 0;    signed short real_part, imaginary_part; /* 定义实部和虚部分量 */    int sqrt_harmonic_sum = 0;  /* 所有谐波分量的平方和 */        if(wave.flag == 1) {        wave.flag = 0;        /* 1. 失能DMA1通道1 以防在中断过程中发生DMA传输 */        DMA_Cmd(ADC_DMA_CH, DISABLE);                /* 2. 将得到的ADC数据处理后放入fft_input_buff数组中 */        for(i=0; i<FFT_POINTS; i++) {            fft_input_buff[i] = ((signed short)(adcValue[i])) << 16;  /* 左移16位 高16位存放实部 */        }        /* 3. 调用FFT库计算 */        cr4_fft_1024_stm32(fft_output_buff, fft_input_buff, FFT_POINTS);                /* 4. 获取信号谐波分量的幅值 */        fft_get_harmonic_data();                /* 5. 获取fft频率数组最大值索引 */        wave.index = findMax_index(&magnitude_buff[1], FFT_POINTS/2-1) + 1;/* 第0个点是FFT直流分量 所以从第一个数据开始取 索引值index需要加1 */                /* 6. 计算波形的频率 幅值 相位等信息 */        wave.freq = wave.sample_freq / (float)FFT_POINTS * (float)(wave.index); /* 频率计算公式: Fs = F0/points*index */        wave.amplitude = magnitude_buff[wave.index] / (FFT_POINTS/2); /* 幅值计算公式: Vpp = buff[index]/(points/2) */                real_part = fft_output_buff[wave.index] >> 16;        /* 计算实部 */        imaginary_part = fft_output_buff[wave.index] & 0xFFFF;/* 计算虚部 */        wave.phase = atan2(real_part, imaginary_part);      /* 相位计算公式: phase = atan2(real, imaginary) */                /* 7. 计算总谐波失真 */        for(i=1; i<FFT_POINTS/2-1; i++) {            sqrt_harmonic_sum += magnitude_buff[i] * magnitude_buff[i];        }        wave.THD = sqrtl(sqrt_harmonic_sum) / magnitude_buff[0];                /* 8. 重新使能DMA1通道1 开始工作 */        DMA_Cmd(ADC_DMA_CH, ENABLE);    }}

故周期运算中增加该部分算法，实现如第7步步骤所示。设计驱动函数中的接口函数如下：

/** * @brief  获取FFT的总谐波失真 单位% * @param  [in] *   @arg  None * @retval FFT总谐波失真 */float fft_get_THD(void){    assert_param(wave.THD);    return (wave.THD);}

后续将逐渐补充其他相关知识点......

# 上一篇：ReaLTaiizor开源.NET winform控件库学

# 下一篇：github加速神器！解决github巨慢的问题，