【STM32F407的DSP教程】第31章 STM32F407實數浮點FFT（支持單精度和雙精度）

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

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

第31章 STM32F407實數浮點FFT（支持單精度和雙精度）

本章主要講解實數浮點FTT，支持單精度和雙精度。

31.1 初學者重要提示

31.2 實數浮點FFT 說明

31.3 單精度函數arm_rfft_fast_f32的使用（含幅頻和相頻）

31.4 雙精度函數arm_rfft_ fast_f64的使用（含幅頻和相頻）

31.5 實驗例程說明(MDK)

31.6 實驗例程說明(IAR)

31.7 總結

31.1 初學者重要提示

與上一章節的復數FFT相比，實數FFT僅需用戶輸入實部即可。輸出結果根據FFT的對稱性，也僅輸出一半的頻譜。

31.2 實數浮點FFT說明

CMSIS DSP庫里面包含一個專門用於計算實數序列的FFT庫，很多情況下，用戶只需要計算實數序列即可。計算同樣點數FFT的實數序列要比計算同樣點數的虛數序列有速度上的優勢。

快速的rfft算法是基於混合基cfft算法實現的。

一個N點的實數序列FFT正變換采用下面的步驟實現：

由上面的框圖可以看出，實數序列的FFT是先計算N/2個實數的CFFT，然后再重塑數據進行處理從而獲得半個FFT頻譜即可（利用了FFT變換后頻譜的對稱性）。

一個N點的實數序列FFT逆變換采用下面的步驟實現：

實數FFT支持浮點，Q31和Q15三種數據類型。

31.3 單精度函數arm_rfft_fast_f32的使用(含幅頻和相頻)

31.3.1 函數說明

函數原型：

void arm_rfft_fast_f32(
  const arm_rfft_fast_instance_f32 * S,
  float32_t * p,
  float32_t * pOut,
  uint8_t ifftFlag)

函數描述：

這個函數用於單精度浮點實數FFT。

函數參數：

1、第1個參數是封裝好的浮點FFT例化，需要用戶先調用函數arm_rfft_fast_init_f32初始化，然后供此函數arm_rfft_fast_f32調用。支持32, 64, 128, 256, 512, 1024, 2048, 4096點FFT。

比如做1024點FFT，代碼如下：

arm_rfft_fast_instance_f32 S;

arm_rfft_fast_init_f32(&S, 1024);

arm_rfft_fast_f32(&S, testInput_f32, testOutput_f32, ifftFlag);

2、第2個參數是實數地址，比如我們要做1024點實數FFT，要保證有1024個緩沖。

3、第3個參數是FFT轉換結果，轉換結果不是實數了，而是復數，按照實部，虛擬，實部，虛部，依次排列。比如做1024點FFT，這里的輸出也會有1024個數據，即512個復位。

4、第4個參數用於設置正變換和逆變換，ifftFlag=0表示正變換，ifftFlag=1表示逆變換。

31.3.2 使用舉例並和Matlab比較

下面通過在開發板上運行這個函數並計算幅頻相應，然后再與Matlab計算的結果做對比。

/*
*********************************************************************************************************
*    函 數 名: arm_rfft_f32_app
*    功能說明: 調用函數arm_rfft_fast_f32計算幅頻和相頻
*    形    參：無
*    返 回 值: 無
*********************************************************************************************************
*/
static void arm_rfft_f32_app(void)
{
    uint16_t i;
    arm_rfft_fast_instance_f32 S;
    
    
    /* 正變換 */
    ifftFlag = 0; 
    
    /* 初始化結構體S中的參數 */
     arm_rfft_fast_init_f32(&S, TEST_LENGTH_SAMPLES);
    
    for(i=0; i<1024; i++)
    {
        /* 波形是由直流分量，50Hz正弦波組成，波形采樣率1024，初始相位60° */
        testInput_f32[i] = 1 + cos(2*3.1415926f*50*i/1024 + 3.1415926f/3);
    }
    
    /* 1024點實序列快速FFT */ 
    arm_rfft_fast_f32(&S, testInput_f32, testOutput_f32, ifftFlag);
    
    /* 為了方便跟函數arm_cfft_f32計算的結果做對比，這里求解了1024組模值，實際函數arm_rfft_fast_f32
       只求解出了512組  
    */ 
     arm_cmplx_mag_f32(testOutput_f32, testOutputMag_f32, TEST_LENGTH_SAMPLES);
    
    
    printf("=========================================\r\n");    
    
    /* 求相頻 */
    PowerPhaseRadians_f32(testOutput_f32, Phase_f32, TEST_LENGTH_SAMPLES, 0.5f);
    
    
    /* 串口打印幅值和相頻 */
    for(i=0; i<TEST_LENGTH_SAMPLES; i++)
    {
        printf("%f, %f\r\n", testOutputMag_f32[i], Phase_f32[i]);
    }
}

運行函數arm_rfft_f32_app可以通過串口打印出計算的模值和相角，下面我們就通過Matlab計算的模值和相角跟arm_rfft_fast_f32計算的做對比。

對比前需要先將串口打印出的數據加載到Matlab中，並給這個數組起名sampledata，加載方法在前面的教程的第13章13.6小結已經講解，這里不做贅述了。Matlab中運行的代碼如下：：

Fs = 1024;               % 采樣率
N  = 1024;               % 采樣點數
n  = 0:N-1;              % 采樣序列
t  = 0:1/Fs:1-1/Fs;      % 時間序列
f = n * Fs / N;          %真實的頻率

%波形是由直流分量，50Hz正弦波正弦波組成
x = 1 + cos(2*pi*50*t + pi/3)   ;  
y = fft(x, N);               %對原始信號做FFT變換
Mag = abs(y);

subplot(2,2,1);
plot(f, Mag); 
title('Matlab計算幅頻響應');
xlabel('頻率');
ylabel('賦值');

subplot(2,2,2);
realvalue = real(y);
imagvalue = imag(y);
plot(f, atan2(imagvalue, realvalue)*180/pi.*(Mag>=200)); 
title('Matlab計算相頻響應');
xlabel('頻率');
ylabel('相角');

subplot(2,2,3);
plot(f, sampledata1);  %繪制STM32計算的幅頻相應
title('STM32計算幅頻響應');
xlabel('頻率');
ylabel('賦值');

subplot(2,2,4);
plot(f, sampledata2);   %繪制STM32計算的相頻相應
title('STM32計算相頻響應');
xlabel('頻率');
ylabel('相角');

運行Matlab后的輸出結果如下：

從上面的對比結果中可以看出，從上面的前512點對比中，我們可以看出兩者的計算結果是相符的Matlab和函數arm_rfft_fast_f32計算的結果基本是一直的。幅頻響應求出的幅值和相頻響應中的求出的初始相角都是沒問題的。

31.4 雙精度函數arm_rfft_fast_f64的使用(含幅頻和相頻)

31.4.1 函數說明

函數原型：

void arm_rfft_fast_f64(
  arm_rfft_fast_instance_f64 * S,
  float64_t * p,
  float64_t * pOut,
  uint8_t ifftFlag)

函數描述：

這個函數用於雙精度浮點實數FFT。

函數參數：

1、第1個參數是封裝好的浮點FFT例化，需要用戶先調用函數arm_rfft_fast_init_f64初始化，然后供此函數arm_rfft_fast_f64調用。支持32, 64, 128, 256, 512, 1024, 2048, 4096點FFT。

比如做1024點FFT，代碼如下：

arm_rfft_fast_instance_f64 S;

arm_rfft_fast_init_f64(&S, 1024);

arm_rfft_fast_f64(&S, testInput_f64, testOutput_f64, ifftFlag);

2、第2個參數是實數地址，比如我們要做1024點實數FFT，要保證有1024個緩沖。

4、第4個參數用於設置正變換和逆變換，ifftFlag=0表示正變換，ifftFlag=1表示逆變換

31.4.2 使用舉例並和Matlab比較

下面通過在開發板上運行這個函數並計算幅頻相應，然后再與Matlab計算的結果做對比。

/*
*********************************************************************************************************
*    函 數 名: arm_rfft_f64_app
*    功能說明: 調用函數arm_rfft_fast_f64計算幅頻和相頻
*    形    參：無
*    返 回 值: 無
*********************************************************************************************************
*/
static void arm_rfft_f64_app(void)
{
    uint16_t i;
    float64_t lX,lY;
    arm_rfft_fast_instance_f64 S;
    
    
    /* 正變換 */
    ifftFlag = 0; 
    
    /* 初始化結構體S中的參數 */
     arm_rfft_fast_init_f64(&S, TEST_LENGTH_SAMPLES);
    
    for(i=0; i<1024; i++)
    {
        /* 波形是由直流分量，50Hz正弦波組成，波形采樣率1024，初始相位60° */
        testInput_f64[i] = 1 + cos(2*3.1415926*50*i/1024 + 3.1415926/3);
    }
    
    /* 1024點實序列快速FFT */ 
    arm_rfft_fast_f64(&S, testInput_f64, testOutput_f64, ifftFlag);
    
    /* 求解模值  */ 
    for (i =0; i < TEST_LENGTH_SAMPLES; i++)
    {
         lX = testOutput_f64[2*i];                    /* 實部*/
        lY = testOutput_f64[2*i+1];                   /* 虛部 */  
        testOutputMag_f64[i] = sqrt(lX*lX+ lY*lY);   /* 求模 */
    }
        
    
    printf("=========================================\r\n");    
    
    /* 求相頻 */
    PowerPhaseRadians_f64(testOutput_f64, Phase_f64, TEST_LENGTH_SAMPLES, 0.5);
    

    /* 串口打印幅值和相頻 */
    for(i=0; i<TEST_LENGTH_SAMPLES; i++)
    {
        printf("%.11f, %.11f\r\n", testOutputMag_f64[i], Phase_f64[i]);
    }    
            
}

運行函數arm_rfft_f64_app可以通過串口打印出計算的模值和相角，下面我們就通過Matlab計算的模值和相角跟arm_rfft_fast_f32計算的做對比。

對比前需要先將串口打印出的數據加載到Matlab中，並給這個數組起名sampledata，加載方法在前面的教程的第13章13.6小結已經講解，這里不做贅述了。Matlab中運行的代碼如下：

Fs = 1024;               % 采樣率
N  = 1024;               % 采樣點數
n  = 0:N-1;              % 采樣序列
t  = 0:1/Fs:1-1/Fs;      % 時間序列
f = n * Fs / N;          %真實的頻率

%波形是由直流分量，50Hz正弦波正弦波組成
x = 1 + cos(2*pi*50*t + pi/3)   ;  
y = fft(x, N);               %對原始信號做FFT變換
Mag = abs(y);

subplot(2,2,1);
plot(f, Mag); 
title('Matlab計算幅頻響應');
xlabel('頻率');
ylabel('賦值');

subplot(2,2,2);
realvalue = real(y);
imagvalue = imag(y);
plot(f, atan2(imagvalue, realvalue)*180/pi.*(Mag>=200)); 
title('Matlab計算相頻響應');
xlabel('頻率');
ylabel('相角');

subplot(2,2,3);
plot(f, sampledata1);  %繪制STM32計算的幅頻相應
title('STM32計算幅頻響應');
xlabel('頻率');
ylabel('賦值');

subplot(2,2,4);
plot(f, sampledata2);   %繪制STM32計算的相頻相應
title('STM32計算相頻響應');
xlabel('頻率');
ylabel('相角');

運行Matlab后的輸出結果如下：

從上面的對比結果中可以看出，從上面的前512點對比中，我們可以看出兩者的計算結果是相符的Matlab和函數arm_rfft_fast_f64計算的結果基本是一直的。幅頻響應求出的幅值和相頻響應中的求出的初始相角都是沒問題的。

31.5 實驗例程說明（MDK）

配套例子：

V5-221_實數浮點FTT(支持單精度和雙精度)

實驗目的：

學習實數浮點FFT，支持單精度浮點和雙精度浮點

實驗內容：

啟動一個自動重裝軟件定時器，每100ms翻轉一次LED2。
按下按鍵K1，串口打印1024點實數單精度FFT的幅頻響應和相頻響應。
按下按鍵K2，串口打印1024點實數雙精度FFT的幅頻響應和相頻響應。

使用AC6注意事項

特別注意附件章節C的問題

上電后串口打印的信息：

波特率 115200，數據位 8，奇偶校驗位無，停止位 1。

RTT方式打印信息：

程序設計：

系統棧大小分配：

硬件外設初始化

硬件外設的初始化是在 bsp.c 文件實現：

/*
*********************************************************************************************************
*    函 數 名: bsp_Init
*    功能說明: 初始化所有的硬件設備。該函數配置CPU寄存器和外設的寄存器並初始化一些全局變量。只需要調用一次
*    形    參：無
*    返 回 值: 無
*********************************************************************************************************
*/
void bsp_Init(void)
{
    /* 
       STM32F407 HAL 庫初始化，此時系統用的還是F407自帶的16MHz，HSI時鍾:
       - 調用函數HAL_InitTick，初始化滴答時鍾中斷1ms。
       - 設置NVIC優先級分組為4。
     */
    HAL_Init();

    /* 
       配置系統時鍾到168MHz
       - 切換使用HSE。
       - 此函數會更新全局變量SystemCoreClock，並重新配置HAL_InitTick。
    */
    SystemClock_Config();

    /* 
       Event Recorder：
       - 可用於代碼執行時間測量，MDK5.25及其以上版本才支持，IAR不支持。
       - 默認不開啟，如果要使能此選項，務必看V5開發板用戶手冊第8章
    */    
#if Enable_EventRecorder == 1  
    /* 初始化EventRecorder並開啟 */
    EventRecorderInitialize(EventRecordAll, 1U);
    EventRecorderStart();
#endif
    
    bsp_InitKey();        /* 按鍵初始化，要放在滴答定時器之前，因為按鈕檢測是通過滴答定時器掃描 */
    bsp_InitTimer();      /* 初始化滴答定時器 */
    bsp_InitUart();    /* 初始化串口 */
    bsp_InitLed();        /* 初始化LED */        
}

主功能：

主程序實現如下操作：

啟動一個自動重裝軟件定時器，每100ms翻轉一次LED2。
按下按鍵K1，串口打印1024點實數單精度FFT的幅頻響應和相頻響應。
按下按鍵K2，串口打印1024點實數雙精度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鍵按下 */
                    arm_rfft_f32_app();
                    break;
                
                case KEY_DOWN_K2:            /* K2鍵按下 */
                    arm_rfft_f64_app();
                    break;
                
                    
                default:
                    /* 其它的鍵值不處理 */
                    break;
            }
        }

    }
}

31.6 實驗例程說明（IAR）