【嵌入式】利用arm-DSP库进行FFT计算，获得信号的频谱、幅值及相位(下)

代码终于要来了，那么咱们就一步一步来。核心代码就三个文件——data.c data.h ffttest.c,其中data.c data.h为数据的准备部分，数据当然来自Matlab。ffttest.c里才是干货。(没看过上的自行翻阅)

一.data.c/data.h数据准备

data.h

#ifndef _DATA_H
#define _DATA_H
#include "arm_math.h"
#include  <stm32f4xx_hal.h>
#define SAMPLE_N    64
#define SAMPLE_FS   3200
#define FFT_DPI     (SAMPLE_FS / SAMPLE_N);

extern float32_t data2handle[SAMPLE_N * 2];
extern float32_t FFT_Output[SAMPLE_N];

extern float32_t xiebo[SAMPLE_N];
extern float32_t phaseA[SAMPLE_N];
extern float32_t phaseB[SAMPLE_N];
extern float32_t phaseC[SAMPLE_N];

#endif

data.c

#include "data.h"

float32_t data2handle[SAMPLE_N * 2] = {0};
float32_t FFT_Output[SAMPLE_N]      = {0};

float32_t xiebo[SAMPLE_N]  = {0,12.0830468476080,23.0625562988315,31.9767233645297,38.1266205974219,41.1592650408757,41.1008388582348,38.3357763183915,33.5355339059327,27.5483591934343,21.2671323190531,15.4955017031488,10.8326010981529,7.59361568157981,5.77782765517063,5.08838180836849,5,4.86346545835348,4.03002514886167,1.97578767574228,-1.59380577304000,-6.67628905966528,-12.9524361960277,-19.8182546598069,-26.4644660940673,-31.9918434767550,-35.5451365280388,-36.4452976726157,-34.2997862737710,-29.0738765919851,-21.1116530786703,-11.1028754443124,-1.40834381901946e-14,11.1028754443124,21.1116530786703,29.0738765919851,34.2997862737710,36.4452976726157,35.5451365280388,31.9918434767550,26.4644660940672,19.8182546598069,12.9524361960277,6.67628905966530,1.59380577303998,-1.97578767574228,-4.03002514886168,-4.86346545835348,-5,-5.08838180836849,-5.77782765517063,-7.59361568157980,-10.8326010981529,-15.4955017031488,-21.2671323190531,-27.5483591934343,-33.5355339059328,-38.3357763183915,-41.1008388582348,-41.1592650408757,-38.1266205974219,-31.9767233645297,-23.0625562988315,-12.0830468476080};

float32_t phaseA[SAMPLE_N] = {0,0.490085701647803,0.975451610080641,1.45142338627231,1.91341716182545,2.35698368412999,2.77785116509801,3.17196642081823,3.53553390593274,3.86505226681368,4.15734806151273,4.40960632174178,4.61939766255643,4.78470167866105,4.90392640201615,4.97592363336099,5,4.97592363336099,4.90392640201615,4.78470167866104,4.61939766255643,4.40960632174178,4.15734806151273,3.86505226681368,3.53553390593274,3.17196642081823,2.77785116509801,2.35698368412999,1.91341716182545,1.45142338627231,0.975451610080641,0.490085701647802,6.12323399573677e-16,-0.490085701647803,-0.975451610080642,-1.45142338627231,-1.91341716182545,-2.35698368412999,-2.77785116509801,-3.17196642081823,-3.53553390593274,-3.86505226681369,-4.15734806151273,-4.40960632174178,-4.61939766255643,-4.78470167866104,-4.90392640201615,-4.97592363336099,-5,-4.97592363336098,-4.90392640201615,-4.78470167866105,-4.61939766255643,-4.40960632174178,-4.15734806151273,-3.86505226681368,-3.53553390593274,-3.17196642081823,-2.77785116509801,-2.35698368412999,-1.91341716182545,-1.45142338627231,-0.975451610080639,-0.490085701647803};
float32_t phaseB[SAMPLE_N] = {-4.33012701892219,-4.55431912460588,-4.73465064747553,-4.86938489638667,-4.95722430686905,-4.99732293738183,-4.98929461619302,-4.93321666042440,-4.82962913144534,-4.67952963378663,-4.48436370766344,-4.24601090763289,-3.96676670145618,-3.64932036348918,-3.29672907550034,-2.91238848433901,-2.50000000000000,-2.06353514902197,-1.60719732651581,-1.13538131517186,-0.652630961100256,-0.163595414108879,0.327015646150717,0.814477366972946,1.29409522551260,1.76125023960617,2.21144345109501,2.64033925325184,3.04380714504361,3.41796151011436,3.75919903739489,4.06423342295808,4.33012701892219,4.55431912460588,4.73465064747553,4.86938489638667,4.95722430686905,4.99732293738183,4.98929461619302,4.93321666042440,4.82962913144534,4.67952963378663,4.48436370766344,4.24601090763289,3.96676670145617,3.64932036348918,3.29672907550034,2.91238848433901,2.50000000000000,2.06353514902197,1.60719732651581,1.13538131517187,0.652630961100256,0.163595414108880,-0.327015646150719,-0.814477366972946,-1.29409522551261,-1.76125023960617,-2.21144345109501,-2.64033925325184,-3.04380714504361,-3.41796151011436,-3.75919903739489,-4.06423342295808};
float32_t phaseC[SAMPLE_N] = {4.33012701892219,4.06423342295808,3.75919903739489,3.41796151011436,3.04380714504360,2.64033925325184,2.21144345109501,1.76125023960617,1.29409522551261,0.814477366972945,0.327015646150716,-0.163595414108881,-0.652630961100257,-1.13538131517187,-1.60719732651581,-2.06353514902197,-2.50000000000000,-2.91238848433901,-3.29672907550034,-3.64932036348918,-3.96676670145617,-4.24601090763289,-4.48436370766344,-4.67952963378663,-4.82962913144534,-4.93321666042440,-4.98929461619302,-4.99732293738183,-4.95722430686905,-4.86938489638667,-4.73465064747553,-4.55431912460588,-4.33012701892220,-4.06423342295808,-3.75919903739489,-3.41796151011436,-3.04380714504360,-2.64033925325184,-2.21144345109501,-1.76125023960617,-1.29409522551260,-0.814477366972945,-0.327015646150718,0.163595414108876,0.652630961100260,1.13538131517187,1.60719732651581,2.06353514902197,2.50000000000000,2.91238848433901,3.29672907550035,3.64932036348918,3.96676670145617,4.24601090763289,4.48436370766344,4.67952963378663,4.82962913144534,4.93321666042440,4.98929461619302,4.99732293738183,4.95722430686905,4.86938489638667,4.73465064747553,4.55431912460588};

无非就是准备了几个数组：
1.float32_t data2handle[SAMPLE_N * 2]，这个用来装处理后的采样数据，所谓处理就是之气那说的虚部填上0，所以对于SAMPLE_N 采样点，这个数组的长度为SAMPLE_N *2;
2.float32_t FFT_Output[SAMPLE_N],这个就是FFT输出的频谱数组，咋看之前都讲过了；
3.float32_t xiebo[SAMPLE_N]，就是那个50Hz、100Hz、200Hz叠加后的信号；
4.float32_t phaseA、B、C[SAMPLE_N]，就是标准三相电的信号；

二.ffttest.c数据处理

ffttest.c

#include "data.h"
#include "arm_const_structs.h"
#include "app_serial.h"
#include "os.h"

void Create_data2handle(float32_t *p)
{
    for(int i = 0;i < SAMPLE_N; i++)
    {
        data2handle[2 * i]     = p[i];
        data2handle[2 * i + 1] = 0;
    }
}

void CalXiebo(void)
{
    
    Create_data2handle(xiebo);
    arm_cfft_f32(&arm_cfft_sR_f32_len64, data2handle, 0, 1);
    arm_cmplx_mag_f32(data2handle, FFT_Output, SAMPLE_N);

    App_SerPrintf("/********************************************************************************/\n"); 
    for(int i = 0,j = 0; i < SAMPLE_N;i++,j++)
    {
        App_SerPrintf("%.8f\t",FFT_Output[i]); 
        if(j >= 5)
        {
            App_SerPrintf("\n"); 
            j = -1;
        }
    }
    App_SerPrintf("\n/********************************************************************************/\n");   
}

void CalPhase(float32_t *phase,char phase_name)
{
    int Index_50Hz = 50 / FFT_DPI;
    float32_t Real,Imaginary,hudu,jiaodu;
    Create_data2handle(phase);
    arm_cfft_f32(&arm_cfft_sR_f32_len64, data2handle, 0, 1);
    arm_cmplx_mag_f32(data2handle, FFT_Output, SAMPLE_N);

    Real      = data2handle[Index_50Hz * 2];
    Imaginary = data2handle[Index_50Hz * 2 + 1];

    hudu = atan2(Imaginary , Real);
    jiaodu = hudu * 180.0 / 3.1415926;
    App_SerPrintf("%c相相位:%f°\n",phase_name,jiaodu);

}

void FFTtest(void)
{
    CalXiebo();
    CalPhase(phaseA,'A');
    CalPhase(phaseB,'B');
    CalPhase(phaseC,'C');
}

啥子，才50来行？？？莫急，这玩意用库就是简单的很。。。函数一个一个来：
1.void Create_data2handle(float32_t *p)；就是产生数组，用于arm_cfft_f32的输入；
2.void CalXiebo(void)；计算谐波，先产生数组，再马上进行两行FFT计算，然后将FFTOutput打印出来，为了美观我是煞费苦心，当然后来发现也不是很美观。。
3.void CalPhase(float32_t *phase,char phase_name)；计算相位。首先获得50Hz的下标，即50/分辨率。然后进行两步计算，其实这里系需要一步就够了，因为计算幅值根本用不到FFTOutput这个数组，第一步arm_cfft_f32已经把傅里叶分解结果存在data2handle里了。然后根据下标找到50Hz分量对应的实部虚部，再进行反正切运算和弧度-角度的转换，出来就是喜闻乐见的角度了；
4.void FFTtest(void)，短小精悍的test函数。。。。

三.万事俱备，只差测试

怎么样。。。在串口助手里打印出来还算能看吧。。
将输出的FFT_Output复制到Matlab里画图
可以看到，在分辨率50Hz的情况下，0代表的是直流分量，1为50Hz分量，2为100Hz分量，4为200Hz分量。之前50、100、200设置的幅值分别是5、30、15，那么用上一篇中所说的计算实际幅值的办法，160/32 = 5；960/32=30；480/32=15,完全OK。

那么再来看相位，A相-90°，B相150°，C相30°，向量图如下图，凑活看吧。跟设置的相位一毛一样。
所以给定一个信号，用FFT就能分解出你想看的分量的全部信息，即频率、幅值、相位，可以说非常给力了。