STM32H7的DSP教程第33章 STM32H7不限制点数FFT实现

Posted 2021-06-30 安富莱电子

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了STM32H7的DSP教程第33章 STM32H7不限制点数FFT实现相关的知识，希望对你有一定的参考价值。

完整版教程下载地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547

第33章 STM32H7不限制点数FFT实现

本章主要讲解不限制点数FFT的实现。

33.1 初学者重要提示

33.2 不限制点数FFT移植

33.3 不限制点数FFT应用说明

33.4 实验例程说明(MDK)

33.5 实验例程说明(IAR)

33.6 总结

33.1 初学者重要提示

由于ARM DSP库限制最大只能4096点，有点不够用，所以整理了个更大点数的。不限制点数，满足2^n即可，n最小值4，即16个点的FFT，而最大值不限。
此FFT算法没有再使用ARM DSP库，重新实现了一个。

33.2 不限制点数FFT移植

33.2.1 移植FFT相关文件

移植下面两个文件fft.c和FFTInc.h到工程：

33.2.2 添加路径

添加路径，大家根据自己的工程来设置即可：

33.3 不限制点数FFT应用说明

33.3.1 支持的点数范围

支持最小16点FFT，最大值不限，但需满足2^n。

33.3.2 函数InitTableFFT

函数原型：

/*
*********************************************************************************************************
*    函 数 名: Int_FFT_TAB
*    功能说明: 正弦和余弦表
*    形    参: 点数
*    返 回 值: 无
*********************************************************************************************************
*/
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384)
float32_t   costab[MAX_FFT_N/2];
float32_t   sintab[MAX_FFT_N/2];
void InitTableFFT(uint32_t n)
{
    uint32_t i;

/* 正常使用下面获取cos和sin值 */
#if 1
    for (i = 0; i < n; i ++ )
    {
         sintab[ i ]=  sin( 2 * PI * i / MAX_FFT_N ); 
        costab[ i ]=  cos( 2 * PI * i / MAX_FFT_N ); 
    }
    
/* 打印出来是方便整理cos值和sin值数组，将其放到内部Flash，从而节省RAM */
#else
    printf("const float32_t sintab[%d] = {\\r\\n", n);
    for (i = 0; i < n; i ++ )
    {
         sintab[ i ]=  sin( 2 * PI * i / MAX_FFT_N ); 
        printf("%.11ff,\\r\\n", sintab[ i ]);
    }    
    printf("};\\r\\n");
    
    printf("const float32_t costab[%d] = {\\r\\n", n);
    for (i = 0; i < n; i ++ )
    {
         sintab[ i ]=  cos( 2 * PI * i / MAX_FFT_N ); 
        printf("%.11ff,\\r\\n", sintab[ i ]);
    }    
    printf("};\\r\\n");
#endif
}
#endif

函数描述：

这个函数用于FFT计算过程中需要用到的正弦和余弦表。对于8192点和16384点已经专门制作了数值表，存到内部Flash，其它点数继续使用的RAM空间，大家可以根据所使用芯片的RAM和Flash大小，选择正弦和余弦值存到RAM还是Flash。

函数参数：

第1个参数是FFT点数。

33.3.3 函数cfft

函数原型：

void cfft(struct compx *_ptr, uint32_t FFT_N )

/*
*********************************************************************************************************
*    函 数 名: cfft
*    功能说明: 对输入的复数组进行快速傅里叶变换（FFT）
*    形    参: *_ptr 复数结构体组的首地址指针struct型 
*             FFT_N 表示点数
*    返 回 值: 无
*********************************************************************************************************
*/
void cfft(struct compx *_ptr, uint32_t FFT_N )
{
    float32_t TempReal1, TempImag1, TempReal2, TempImag2;
    uint32_t k,i,j,z;
    uint32_t Butterfly_NoPerColumn;                    /* 每级蝶形的蝶形组数 */
    uint32_t Butterfly_NoOfGroup;                    /* 蝶形组的第几组 */
    uint32_t Butterfly_NoPerGroup;                    /* 蝶形组的第几个蝶形 */
    uint32_t ButterflyIndex1,ButterflyIndex2,P,J;
    uint32_t L;
    uint32_t M;

    z=FFT_N/2;    /* 变址运算，即把自然顺序变成倒位序，采用雷德算法 */
    
    for(i=0,j=0;i<FFT_N-1;i++)        
    {
        /* 
          如果i<j,即进行变址 i=j说明是它本身，i>j说明前面已经变换过了，不许再变化，注意这里一般是实数
 有虚数部分 设置成结合体 
        */
        if(i<j)                                        
        {                                            
            TempReal1  = _ptr[j].real;               
            _ptr[j].real= _ptr[i].real;
            _ptr[i].real= TempReal1;
        }
         
        k=z;         /*求j的下一个倒位序 */
        
        while(k<=j)  /* 如果k<=j,表示j的最高位为1 */  
        {           
            j=j-k; /* 把最高位变成0 */
            k=k/2; /* k/2，比较次高位，依次类推，逐个比较，直到某个位为0，通过下面那句j=j+k使其变为1 */
        }
        
        j=j+k;     /* 求下一个反序号，如果是0，则把0改为1 */
    }
    
/* 第L级蝶形(M)第Butterfly_NoOfGroup组(Butterfly_NoPerColumn)第J个蝶形(Butterfly_NoPerGroup)****** */
/* 蝶形的组数以2的倍数递减Butterfly_NoPerColumn，每组中蝶形的个数以2的倍数递增Butterfly_NoPerGroup */
/* 在计算蝶形时，每L列的蝶形组数,一共有M列，每组蝶形中蝶形的个数,蝶形的阶数(0,1,2.....M-1) */
    Butterfly_NoPerColumn = FFT_N;                             
    Butterfly_NoPerGroup = 1;    
    M = log2(FFT_N);
    for ( L = 0;L < M; L++ )                                 
    {
        Butterfly_NoPerColumn >>= 1;        /* 蝶形组数 假如N=8，则(4,2,1) */
        
        //第L级蝶形 第Butterfly_NoOfGroup组    （0,1，....Butterfly_NoOfGroup-1）                    
        for ( Butterfly_NoOfGroup = 0;Butterfly_NoOfGroup < Butterfly_NoPerColumn;Butterfly_NoOfGroup++ )
        {  
            /* 第 Butterfly_NoOfGroup 组中的第J个 */
            for ( J = 0;J < Butterfly_NoPerGroup;J ++ )        
            {   /* 第 ButterflyIndex1 和第 ButterflyIndex2 个元素作蝶形运算,WNC */
/* (0,2,4,6)(0,1,4,5)(0,1,2,3) */
/* 两个要做蝶形运算的数相距Butterfly_NoPerGroup，ge=1,2,4 */                        
/* 这里相当于P=J*2^(M-L),做了一个换算下标都是N (0,0,0,0)(0,2,0,2)(0,1,2,3) */
                ButterflyIndex1 = ( ( Butterfly_NoOfGroup * Butterfly_NoPerGroup ) << 1 ) + J;
                ButterflyIndex2 = ButterflyIndex1 + Butterfly_NoPerGroup;
                P = J * Butterfly_NoPerColumn;                
                
                /* 计算和转换因子乘积 */
                TempReal2 = _ptr[ButterflyIndex2].real * costab[ P ] +  _ptr[ButterflyIndex2].imag * 
sintab[ P ];
                TempImag2 = _ptr[ButterflyIndex2].imag * costab[ P ] -  _ptr[ButterflyIndex2].real * 
sintab[ P ] ;
                TempReal1 = _ptr[ButterflyIndex1].real;
                TempImag1 = _ptr[ButterflyIndex1].imag;
                
                /* 蝶形运算 */
                _ptr[ButterflyIndex1].real = TempReal1 + TempReal2;    
                _ptr[ButterflyIndex1].imag = TempImag1 + TempImag2;
                _ptr[ButterflyIndex2].real = TempReal1 - TempReal2;
                _ptr[ButterflyIndex2].imag = TempImag1 - TempImag2;
            }
        } 
        
        Butterfly_NoPerGroup<<=1; /* 一组中蝶形的个数(1,2,4) */
    }
}

函数描述：

这个函数用于复数FFT变换。

函数参数：

第1个参数是复数格式。
第2个参数是FFT点数，最小值16，最大值不限，满足满足2^n即可。

33.3.4 FFT幅频响应举例

下面通过函数cff将正弦波做16384点FFT变换：

/*
*********************************************************************************************************
*    函 数 名: cfft_f32_mag
*    功能说明: 计算幅频
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
static void cfft_f32_mag(void)
{
    uint32_t i;
    
    /* 计算一批sin，cos系数 */
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384)
    InitTableFFT(MAX_FFT_N);
#endif
    
    for(i=0; i<MAX_FFT_N; i++)
    {
        /* 波形是由直流分量，500Hz正弦波组成，波形采样率MAX_FFT_N，初始相位60° */
        s[i].real = 1 + cos(2*3.1415926f*500*i/MAX_FFT_N + 3.1415926f/3);        
        s[i].imag = 0;
    }
    
    /* MAX_FFT_N点快速FFT */ 
    cfft(s, MAX_FFT_N);
    
    /* 计算幅频 */ 
    for(i=0; i<MAX_FFT_N; i++)
    {
        arm_sqrt_f32((float32_t)(s[i].real *s[i].real+ s[i].imag*s[i].imag ), &s[i].real); 
    }
    
    /* 串口打印求解的幅频 */
    for(i=0; i<MAX_FFT_N; i++)
    {
        printf("%f\\r\\n", s[i].real);
    }
}

运行函数cfft_f32_mag可以通过串口打印FFT结果：

从上面的结果中可以出直流分量和正弦波幅值都是没有问题的。

33.3.5 FFT相频响应举例

下面通过函数cff将正弦波做16384点FFT变换：

/*
*********************************************************************************************************
*    函 数 名: PowerPhaseRadians_f32
*    功能说明: 求相位
*    形    参：_usFFTPoints  复数个数，每个复数是两个float32_t数值
*             _uiCmpValue  比较值，需要求出相位的数值
*    返 回 值: 无
*********************************************************************************************************
*/
void PowerPhaseRadians_f32(uint16_t _usFFTPoints, float32_t _uiCmpValue)        
{
    float32_t lX, lY;
    uint32_t i;
    float32_t phase;
    float32_t mag;
    
    
    for (i=0; i <_usFFTPoints; i++)
    {
        lX= s[i].real;       /* 实部 */
        lY= s[i].imag;    /* 虚部 */ 
        
         phase = atan2f(lY, lX);                            /* atan2求解的结果范围是(-pi, pi], 弧度制 */
        arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag);   /* 求模 */
        
        if(_uiCmpValue > mag)
        {
            s[i].real = 0;            
        }
        else
        {
            s[i].real= phase* 180.0f/3.1415926f;   /* 将求解的结果由弧度转换为角度 */
        }
    }
}

/*
*********************************************************************************************************
*    函 数 名: cfft_f32_phase
*    功能说明: 计算相频
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
static void cfft_f32_phase(void)
{
    uint32_t i;
    

    /* 计算一批sin，cos系数 */
#if (MAX_FFT_N != 8192) && (MAX_FFT_N != 16384)
    InitTableFFT(MAX_FFT_N);
#endif
    
    for(i=0; i<MAX_FFT_N; i++)
    {
        /* 波形是由直流分量，500Hz正弦波组成，波形采样率MAX_FFT_N，初始相位60° */
        s[i].real = 1 + cos(2*3.1415926f*500*i/MAX_FFT_N + 3.1415926f/3);        
        s[i].imag = 0;
    }
    
    /* MAX_FFT_N点快速FFT */ 
    cfft(s, MAX_FFT_N);
    
    /* 求相频 */
    PowerPhaseRadians_f32(MAX_FFT_N, 0.5f);
    
    /* 串口打印求解相频 */
    for(i=0; i<MAX_FFT_N; i++)
    {
        printf("%f\\r\\n", s[i].real);
    }
    
}

运行函数cfft_f32_phase可以通过串口打印FFT结果：

从上面的结果中可以出计算的初始相位是没有问题的。

33.4 实验例程说明（MDK）

配套例子：

V7-223_不限制点数FFT实现

实验目的：

学习不限制点数FFT。

实验内容：

启动一个自动重装软件定时器，每100ms翻转一次LED2。
按下按键K1，串口打印16384点FFT的幅频响应。
按下按键K2，串口打印16384点FFT的相频响应。

使用AC6注意事项

特别注意附件章节C的问题

上电后串口打印的信息：

波特率 115200，数据位 8，奇偶校验位无，停止位 1。

RTT方式打印信息：

程序设计：

系统栈大小分配：

RAM空间用的DTCM：

硬件外设初始化

硬件外设的初始化是在 bsp.c 文件实现：

/*
*********************************************************************************************************
*    函 数 名: bsp_Init
*    功能说明: 初始化所有的硬件设备。该函数配置CPU寄存器和外设的寄存器并初始化一些全局变量。只需要调用一次
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
void bsp_Init(void)
{
    /* 配置MPU */
    MPU_Config();
    
    /* 使能L1 Cache */
    CPU_CACHE_Enable();

    /* 
       STM32H7xx HAL 库初始化，此时系统用的还是H7自带的64MHz，HSI时钟:
       - 调用函数HAL_InitTick，初始化滴答时钟中断1ms。
       - 设置NVIC优先级分组为4。
     */
    HAL_Init();

    /* 
       配置系统时钟到400MHz
       - 切换使用HSE。
       - 此函数会更新全局变量SystemCoreClock，并重新配置HAL_InitTick。
    */
    SystemClock_Config();

    /* 
       Event Recorder：
       - 可用于代码执行时间测量，MDK5.25及其以上版本才支持，IAR不支持。
       - 默认不开启，如果要使能此选项，务必看V7开发板用户手册第8章
    */    
#if Enable_EventRecorder == 1  
    /* 初始化EventRecorder并开启 */
    EventRecorderInitialize(EventRecordAll, 1U);
    EventRecorderStart();
#endif
    
    bsp_InitKey();        /* 按键初始化，要放在滴答定时器之前，因为按钮检测是通过滴答定时器扫描 */
    bsp_InitTimer();      /* 初始化滴答定时器 */
    bsp_InitUart();    /* 初始化串口 */
    bsp_InitExtIO();    /* 初始化FMC总线74HC574扩展IO. 必须在 bsp_InitLed()前执行 */    
    bsp_InitLed();        /* 初始化LED */    
}

MPU配置和Cache配置：

数据Cache和指令Cache都开启。配置了AXI SRAM区（本例子未用到AXI SRAM），FMC的扩展IO区。

/*
*********************************************************************************************************
*    函 数 名: MPU_Config
*    功能说明: 配置MPU
*    形    参: 无
*    返 回 值: 无
*********************************************************************************************************
*/
static void MPU_Config( void )
{
    MPU_Region_InitTypeDef MPU_InitStruct;

    /* 禁止 MPU */
    HAL_MPU_Disable();

    /* 配置AXI SRAM的MPU属性为关闭读Cache和写Cache */
    MPU_InitStruct.Enable           = MPU_REGION_ENABLE;
    MPU_InitStruct.BaseAddress      = 0x24000000;
    MPU_InitStruct.Size             = MPU_REGION_SIZE_512KB;
    MPU_InitStruct.AccessPermission = MPU_REGION_FULL_ACCESS;
    MPU_InitStruct.IsBufferable     = MPU_ACCESS_NOT _BUFFERABLE;
    MPU_InitStruct.IsCacheable      = MPU_ACCESS_NOT _CACHEABLE;
    MPU_InitStruct.IsShareable      = MPU_ACCESS_NOT_SHAREABLE;
    MPU_InitStruct.Number           = MPU_REGION_NUMBER0;
    MPU_InitStruct.TypeExtField     = MPU_TEX_LEVEL0;
    MPU_InitStruct.SubRegionDisable = 0x00;
    MPU_InitStruct.DisableExec      = MPU_INSTRUCTION_ACCESS_ENABLE;

    HAL_MPU_ConfigRegion(&MPU_InitStruct);
    
    
    /* 配置FMC扩展IO的MPU属性为Device或者Strongly Ordered */
    MPU_InitStruct.Enable           = MPU_REGION_ENABLE;
    MPU_InitStruct.BaseAddress      = 0x60000000;
    MPU_InitStruct.Size             = ARM_MPU_REGION_SIZE_64KB;    
    MPU_InitStruct.AccessPermission = MPU_REGION_FULL_ACCESS;
    MPU_InitStruct.IsBufferable     = MPU_ACCESS_BUFFERABLE;
    MPU_InitStruct.IsCacheable      = MPU_ACCESS_NOT_CACHEABLE;    
    MPU_InitStruct.IsShareable      = MPU_ACCESS_NOT_SHAREABLE;
    MPU_InitStruct.Number           = MPU_REGION_NUMBER1;
    MPU_InitStruct.TypeExtField     = MPU_TEX_LEVEL0;
    MPU_InitStruct.SubRegionDisable = 0x00;
    MPU_InitStruct.DisableExec      = MPU_INSTRUCTION_ACCESS_ENABLE;
    
    HAL_MPU_ConfigRegion(&MPU_InitStruct);

    /*使能 MPU */
    HAL_MPU_Enable(MPU_PRIVILEGED_DEFAULT);
}

/*
*********************************************************************************************************
*    函 数 名: CPU_CACHE_Enable
*    功能说明: 使能L1 Cache
*    形    参: 无
*    返 回 值: 无
*********************************************************************************************************
*/
static void CPU_CACHE_Enable(void)
{
    /* 使能 I-Cache */
    SCB_EnableICache();

    /* 使能 D-Cache */
    SCB_EnableDCache();
}

主功能：

主程序实现如下操作：

启动一个自动重装软件定时器，每100ms翻转一次LED2。
按下按键K1，串口打印16384点FFT的幅频响应。
按下按键K2，串口打印16384点FFT的相频响应。

/*
*********************************************************************************************************
*    函 数 名: main
*    功能说明: c程序入口
*    形    参: 无
*    返 回 值: 错误代码(无需处理)
*********************************************************************************************************
*/
int main(void)
{
    uint8_t ucKeyCode;        /* 按键代码 */
    

    bsp_Init();        /* 硬件初始化 */
    PrintfLogo();    /* 打印例程信息到串口1 */

    PrintfHelp();    /* 打印操作提示信息 */
    

    bsp_StartAutoTimer(0, 100);    /* 启动1个100ms的自动重装的定时器 */

    /* 进入主程序循环体 */
    while (1)
    {
        bsp_Idle();        /* 这个函数在bsp.c文件。用户可以修改这个函数实现CPU休眠和喂狗 */
        

        if (bsp_CheckTimer(0))    /* 判断定时器超时时间 */
        {
            /* 每隔100ms 进来一次 */
            bsp_LedToggle(4);    /* 翻转LED2的状态 */   
        }
        
        ucKeyCode = bsp_GetKey();    /* 读取键值, 无键按下时返回 KEY_NONE = 0 */
        if (ucKeyCode != KEY_NONE)
        {
            switch (ucKeyCode)
            {
                case KEY_DOWN_K1:            /* K1键按下 */
                    cfft_f32_mag();
                    break;
                
                case KEY_DOWN_K2:            /* K2键按下 */
                    cfft_f32_phase();
                    break;
                
                    
                default:
                    /* 其它的键值不处理 */
                    break;
            }
        }

    }
}

33.5 实验例程说明（IAR）