电力系统中往往掺杂谐波,而FFT可以将谐波检测出来,具有较大的实用价值。今天主要讲一下在STM32中如何利用dsp库进行快速傅里叶计算,从而得出信号的频谱幅值以及相位。

一、Matlab简单搭建

1.谐波检测搭建

通过50Hz正弦波叠加100Hz、200Hz的正弦波得到叠加后畸变的波形,这个波形主要用来进行FFT运算进行谐波检测。模型很简单如图1所示,三个正弦波相位相同、频率不同,幅值分别为5、15、30。图2中上边是叠加后的畸变的波形,下边是三个信号放在一张图时的波形。
图1图2

2.相位分析搭建

这个很简单,就是搞三个频率、幅值相同,相位分别差120°的正弦波。其中以A的相位为基准,B滞后于A120°,C滞后于B120°。放一起的波形如图2所示,就是电网中典型的三相电。
在这里插入图片描述在这里插入图片描述

3.Matlab数据导入单片机中

这里直接采用在单片机中定义数组,来存储Matlab波形信号。还要提一点,在simulink中设置的仿真步长为0.0003125,其对应的采样值为3200Hz,采样时间设置为0.02秒,即50Hz正弦波的一个周期。那么在这一个周期内,共由64个采样点组成。
将workspace里的数据复制到word中,然后全选,用“,”来替换换行符“^p”,就得到了可以直接在C语言里使用的数组,如下图。
在这里插入图片描述

二、FFT软件编程

1.dsp库添加

如果你使用的是STM32F4,那么巧了,直接从我网盘里下载就好了
链接:https://pan.baidu.com/s/1W5GhIk0GlLFG2clPl-vo0g
提取码:i5mm
如果不是,请自行去ST官网下载。。。
针对FFT运算,只需要把en.stm32f4_dsp_stdperiph_lib\STM32F4xx_DSP_StdPeriph_Lib_V1.8.0\Libraries\CMSIS中的以下7个文件添加到工程中就好了。再就是根据自己的芯片添加宏定义,比如我是M4,我就需要在预编译中加入ARM_MATH_CM4宏定义。在这里插入图片描述

---------2020.7更新,dsp库的添加方法-------------

时隔一年,后来项目中又要用到dsp库,发现并不用像上面提到的,把多个源码添加到工程文件里。只需要将dsp的lib库添加到工程里,然后在需要调用库函数的c文件包含arm_math.h即可。lib文件在CMSIS\Lib\ARM文件夹下,一般单片机是小端模式,所以添加arm_cortexM4lf_math.lib即可。同时别忘了添加ARM_MATH_CM4宏定义。

2.FFT两个重要的接口

1.arm_cfft_f32(const arm_cfft_instance_f32 * S,float32_t * p1,uint8_t ifftFlag,uint8_t bitReverseFlag);这个就是快速傅里叶变换的主要接口,第一个参数可以理解为你输入到FFT里的采样点的个数;第二个参数为输入数组;第三个参数为正反变换,一般使用填0;第四个参数为位反转使能,一般使用填1。对于后两个参数,因为我们的采样数组是原始的采样值,所以一般都需要进行位反转,因为FFT计算要求就需要将输入值进行反转以适应算法。经过这个函数后,我们输入的数组就被傅里叶分解了,数组中每两个元素代表一个数,第一个元素为实部、第二个元素为虚部。
2.arm_cmplx_mag_f32( float32_t * pSrc,float32_t * pDst,uint32_t numSamples);这个是输出频谱的函数。第一个参数为上一个函数傅里叶分解后的数组,第二个参数为频谱的输出数组,最后是采样点的个数。
3.举例:

arm_cfft_f32(&arm_cfft_sR_f32_len64, data2handle, 0, 1);
arm_cmplx_mag_f32(data2handle, FFT_Output, SAMPLE_N);

这里我采用的是64点采样,data2handle是处理后的输入数组,FFT_Output为频谱输出数组,SAMPLE_N为采样点个数,即64。为啥说data2handle是处理后的输入数组呢,因为arm_cfft_f32的输入数组要求为a+bi的形式,即有实部也有虚部,每两个元素组成一个数,下表为偶数的元素为实部,奇数的为虚部。而我们单片机,或者matlab中导出的数据肯定都是实数,故在带入arm_cfft_f32前,需要将数组进行一点小小的处理,让其虚部全为0即可,代码如下,传入函数的参数即为AD采样或者matab导出的数组。

void Create_data2handle(float32_t *p)
{
    for(int i = 0;i < SAMPLE_N; i++)
    {
        data2handle[2 * i]     = p[i];
        data2handle[2 * i + 1] = 0;
    }
}
3.FFT输出的频谱数组咋看

输出的频谱横轴为频率,纵轴为幅度(注意,这个幅度不是幅值,得到实际幅值还需要进一步计算)。所以搞明白频率和幅度咋看就OK了。前边我们将频谱输出到了FFTOutput中了,那咋看呢?
1.频率:首先要明确一个概念->分辨率,它的定义就是频谱中每两个相邻频率之间的频率间隔。计算公式为分辨率=采样频率/采样点数,比如我提到过我的采样频率为3200Hz,采样了一个周波 64个点,那么我的分辨率就是3200/64=50Hz,这就意味着在FFTOutput这个数组中,下标0对应的元素就是0Hz(也就是直流分量)的幅度,下标1对应的就是50Hz的幅度,下标2对应100Hz的幅度,以此类推。如果采样点为3200,那么我们就可以看0、1、2、3…Hz的频率了。
2.幅度:刚刚说了幅度就是输出数组中的元素,但是这并不是实际幅值,如何得到实际幅值呢。除了直流分量要除以N,其它频率的幅度除以N/2就能得到各频率分量的实际幅值了,别忘了N是采样点数哈(别问我为啥,站在巨人的肩膀上,辅之实践得真知)。。。

4.FFT分解后的相位咋看

先复习一个大学工程数学学习的知识。。。对于一个数a+bi,其模值为根号下a方加b方,其相位(即a+bi这个向量与实轴的夹角)为acrtan(b/a)。好了,之前不是说了经过arm_cfft_f32函数后就能得到FFT分解后的结果了吗,还说了两个元素代表一个数,第一个元素为实部,第二个元素为虚部。那么问题就迎刃而解了。比如你要看50Hz分量的相位,那么就去找50Hz对应的那俩元素,直接acrtan(b/a)得出的就是弧度。如果你觉得弧度不好看,那么你就弧度*180/pi,得出的就是我们喜闻乐见的角度了。

写太多了。。分两部分吧。这部分主要讲了基本的理论知识,下一部分主要是代码实现及验证。

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐