stm32f1单片机上用FFT测量信号频率（高精度、过程详细）_stm32 fft

在STM32F1系列单片机上面实现FFT

最近需要做一台基于stm32的示波器。如果测量信号参数，用单片机上面一些传统的方法局限性还是比较大，就开始琢磨直接上FFT。本文将以一个实例来介绍如何使用STM32提供的DSP库函数进行FFT。

开始是下载stm32的DSP库，提供一个下载地址：
链接：https://pan.baidu.com/s/1M28W7WXNsO3QVfkZunh7rg
提取码：s7kb
下载并解压文件：

可以看到官方提供的是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点数
1
2
3

在信号与系统或数字信号处理课程里面经常遇到FFT相关运算题目，根据经验，时域信号通常是实数，但是计算后得到的频域数据往往是虚数，同样按照官方库的说明，输入、输出两个数组都应该是32位的，高16位存储实部，低16位存储虚部。因此在定义数组时：

//   long是32位
long InBufArray[NPT]={0};         //定义输入数组
long OutBufArray[NPT/2];          //定义输出数组
long MagBufArray[NPT/2];          //幅值
1
2
3
4

此库只能单纯的进行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[]为计算输出的幅值数组
    }
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

有了上面的铺垫后就可以用单片机进行数据采集，然后用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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

注意定时器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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

此处配置了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);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

由理论知识可知，幅值最大值的横坐标对应信号频率，纵坐标对应幅度。假设求得幅度值最大的数组为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%。

分析造成误差原因：

1. FFT的分辨率，前面讲到FFT分辨率=Fs/N，要减小误差可增大N，实验用到的N=1024，Fs越大分辨率就越低，误差就越大，因此在测量低频信号时应该降低采样频率，最低要大于信号频率的两倍。
2. ADC读取数据会耗时，定时器触发AD转换时，ADC并非瞬间获取数据，而是需要有短暂的采样时间和转换时间，这样采样频率就不是理想状态下的采样频率，因此可能造成误差。

FFT的优缺点：

1. 测量信号精度高，且最高可以测1M频率的信号，这是我上一篇文章所讲方法《用ADC和定时器测信号频率》无法比拟的，但信号频率很高，采样频率很高时，主循环将受到很大影响，因为DMA和系统总线是复用的，DMA最大能占到二分之一带宽，至少要留给系统二分之一带宽。当采样频率到2M后，由于定时器触发太频繁，导致DMA搬运ADC数据使用太多带宽，会出现主循环进不去的问题。
2. 可以测量很微弱的信号，上一篇文章中讲到：信号幅度没有2V，单片机是无法检测到上升沿的，2V以下信号没办法测量频率，而FFT不一样，最低可以检测幅值为20mV的信号，包括正弦波、方波、锯齿波等等。
3. 测量低频信号很耗时。

关于幅度部分的实验还需等到有函数发生器再进行!!!

(水平有限，讲的不对的地方还望见谅并欢迎指出！)

——————————————————————————————————————————————————
——————————————————————————————————————————————————
代码获取链接：FFT测频率代码

一些收到源码的同学怀疑代码有问题，在这里申明一下：此程序就有很多人收到后成功测量出信号的频率，并且这个只能用在stm32F1系列单片机上，需要适当的调试才能达到测试效果。