【STM32H7的DSP教程】第33章 STM32H7不限制點數FFT實現

本文轉載自查看原文 2021-06-30 17:02 145 第2版DSP和CMSIS-NN神經網絡教程

完整版教程下載地址：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）