stm32f1单片机上用FFT测量信号频率(高精度、过程详细)

时间: 2023-10-04 admin IT培训

stm32f1单片机上用FFT测量信号频率(高精度、过程详细)

stm32f1单片机上用FFT测量信号频率(高精度、过程详细)

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

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

开始是下载stm32的DSP库,提供一个下载地址:
链接:
提取码: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点数

在信号与系统或数字信号处理课程里面经常遇到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);elseMagBufArray[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模式,用来触发ADC1TIM_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%。

分析造成误差原因:

  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. 测量低频信号很耗时。

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

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

——————————————————————————————————————————————————
——————————————————————————————————————————————————
2021年1月1日更新(中间一年时间因为考研去了,很抱歉很多私信和评论都没有回复)

源代码永久链接: .
提取码:aqd2

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