基于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 \r\n location = %d \r\n",maxvalue,Index);
printf("Fmaxvalue = %f \r\n Amplitude = %f \r\n DC_Component = %f \r\n Virtual_value = %f \r\n Res = %f \r\n ",Freq,Amplitude,DC_Component,Virtual_value,res);
至于有效值dsp库中也有支持的函数可以计算，但是我暂时没用到，所以这次就先不写的，挖个坑留到下次写
具体的工程文件我放在 github的仓库了，需要的同学可以自行取用，觉得有帮助的话帮我点个赞吧~
如果发现错误也希望在评论区指出，谢谢大家
1526
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
