FPU DSP软件库（FPU DSP Software Library）详解快速傅里叶变换

德州仪器TMS320C28x浮点单元数字信号处理（FPU DSP）库是一个优化的信号处理例程的集合，为C2000设备编写，支持单精度浮点单元（FPU32），带有三角测量数学单元FPU32（TMU类型0），或双精度FPU（FPU64）。
这些功能使C/C++程序员能够充分利用上述硬件加速器来加快计算时间。FPU DSP Software Library usres guide提供了对库中包含的每个功能的描述。
其中5.4节介绍了实现复数FFT的模块

本文介绍用CFFT_f32()和CFFT_f32t()实现FFT的方法

1.数据结构

官方定义结构体CFFT_F32_STRUCT

结构体成员如下

IInPtr为指向输入缓冲区的指针
OutPtr为指向输出缓冲区的指针

CoefPtr为指向转化因子数组的指针

CurrentInPtr是指向每一阶输入数组的指针
CurrentOutPtr是指向每一阶输出数组的指针。

stage为傅里叶变换的阶数
FFTSize为FFT变换的长度,且有FFTSize=2^stage

2.函数void CFFT_f32 (CFFT_F32_STRUCT_Handle hndCFFT_F32)

该例程以乒乓的方式使用两个缓冲区，即在每个FFT阶段之后，输出和输入缓冲区分别成为下一阶段的输入和输出缓冲区。使用两个指针，CurrentInPtr和CurrentOutPtr来跟踪切换。用户可以通过查看CurrentOutPtr来确定最终输出的地址。
注意： 该例程要求使用两个缓冲区，每个缓冲区的大小为2N（32位float）以进行计算；输入缓冲区必须与4N个字（16位）的内存地址对齐。[^1]如果无法对齐，用户可以使用替代函数CFFT_f32u，尽管速度较慢。
这个函数不能重入，因为它使用全局变量来存储某些参数。

3.void CFFT_f32t (CFFT_F32_STRUCT_Handle hndCFFT_F32)

这个例程计算一个N-pt（N = 2^n，n = 5： 10）复数输入的32位浮点FFT。这个函数将在继续使用FFT之前以位反转的格式重新排序输入。该例程以乒乓的方式使用两个缓冲区，即在每个FFT阶段之后，输出和输入缓冲区分别成为下一阶段的输入和输出缓冲区。使用两个指针，CurrentInPtr和CurrentOutPtr来跟踪切换。用户可以通过查看CurrentOutPtr来确定最终输出的地址。

void CFFT_f32()和void CFFT_f32t()的区别是void CFFT_f32t()使用一个预先生成的旋转因子表。

注意：以上2和3的描述来自于FPU DSP Software Library usres guide。在运行例程的过程中发现用户手册与例程中的描述有出入。

4.例程

中文是我自己加的注释

#define CFFT_STAGES     7//7阶傅里叶变换
#define CFFT_SIZE       (1 << CFFT_STAGES) 
#define EPSILON         0.01

//*****************************************************************************
// globals
//*****************************************************************************
#ifdef __cplusplus
#pragma DATA_SECTION("CFFTdata1")//c++
#else
#pragma DATA_SECTION(CFFTin1Buff,"CFFTdata1")//c语言
#endif //__cplusplus
//! \brief FFT Calculation Buffer
//! \note The input buffer needs to be aligned to a 4N word boundary
//! \note If the number of FFT stages is odd, the result of the FFT will
//! be written to this buffer
//! 
float CFFTin1Buff[CFFT_SIZE*2];//CFFTin1Buff是作为输入缓冲区//需要进行4字节对齐的操作//如果FFT的阶数是奇数，那么FFT的最终结果会写入此缓冲区

#ifdef __cplusplus
#pragma DATA_SECTION("CFFTdata2")
#else
#pragma DATA_SECTION(CFFTin2Buff,"CFFTdata2")
#endif //__cplusplus
//! \brief Magnitude Calculation Buffer//
//!
float CFFTin2Buff[CFFT_SIZE*2];               

#ifdef __cplusplus
#pragma DATA_SECTION("CFFTdata3")
#else
#pragma DATA_SECTION(CFFToutBuff,"CFFTdata3")
#endif //__cplusplus
//! \brief FFT Calculation Buffer
//! \note If the number of FFT stages is even, the result of the FFT will
//! be written to this buffer//
//!
float CFFToutBuff[CFFT_SIZE*2];//CFFToutBuff是作为输出缓冲区//如果FFT的阶数是偶数，那么FFT的最终结果会写入此缓冲区

#ifdef __cplusplus
#pragma DATA_SECTION("CFFTdata4")
#else
#pragma DATA_SECTION(CFFTF32Coef,"CFFTdata4")
#endif //__cplusplus
//! \brief Twiddle Factors
//!
float CFFTF32Coef[CFFT_SIZE];//进行FFT计算所需要的旋转因子

float CFFTgoldenOut[CFFT_SIZE*2] = {//matlab计算fft的结果，用来与DSP的结果进行验证
    #include "data_output_1.h"
};

float CFFTgoldenMagnitude[CFFT_SIZE] = {//matlab计算fft的幅值结果，用来与DSP的结果进行验证
    #include "data_output_2.h"
};

float CFFTgoldenPhase[CFFT_SIZE] = {//matlab计算fft的相位结果，用来与DSP的结果进行验证
    #include "data_output_3.h"
};

float RadStep = 0.1963495408494f;
float Rad = 0.0f;

//! \brief Object of the structure CFFT_F32_STRUCT
//!
CFFT_F32_STRUCT cfft;//官方定义结构体

//! \brief Handle to the CFFT_F32_STRUCT object
//!
CFFT_F32_STRUCT_Handle hnd_cfft = &cfft;

uint16_t pass = 0;
uint16_t fail = 0;

#ifdef USE_TABLES//如果用到预先生成的旋转因子表，需要进行下列操作
//Linker defined variables
extern uint16_t  FFTTwiddlesRunStart;
extern uint16_t  FFTTwiddlesLoadStart;
extern uint16_t  FFTTwiddlesLoadSize;
#endif //USE_TABLES
//*****************************************************************************
// Function Prototypes
//*****************************************************************************

//*****************************************************************************
// function definitions
//*****************************************************************************
//!
//! \brief main routine for the Complex FFT example
//! \return returns a 1
//!
//! The test functions include CFFT_f32,CFFT_f32_mag and CFFT_f32_phase.
//! Data section alignment (#pragma ...) is not necessary for CFFT_f32u 
//! but necessary for testing CFFT_f32.//在CMD文件中实现地址对齐
//! 
//! Minimum CFFT_Stages is 3. When CFFT_Stages is more than 9, the
//! quantization error would be significant. The results can be compared 
//! against MATLAB code under 2837x_CFFT/matlab/CFFTforC28xNew.
//!
//! CFFT_F32_STRUCT is a structure defined as:
//!
//!     typedef struct {
//!        float    *InPtr;
//!        float    *OutPtr;
//!        float    *CoefPtr;
//!        float    *CurrentInPtr;
//!        float    *CurrentOutPtr;
//!        short    Stages;
//!        uint16_t  FFTSize;
//!       }CFFT_F32_STRUCT;
//!   
//!  Watch Variables:
//!
//! -# InPtr           Input/output or middle stage of ping-pong buffer
//! -# OutPtr          Output or middle stage of ping-pong buffer
//! -# CurrentInPtr    Output buffer for CFFT result
//! -# CurrentOutPtr   N-1 stage CFFT result/Magnitude/Phase output buffer
//! -# CoefPtr         Twiddle factor buffer
//!
//! ------------------------------------------------------------------------
//!             | Stage 3     | Stage 4     |..| StageN(odd) | StageN(even)
//! ------------------------------------------------------------------------
//! InPtr(buf1) |CurrentInPtr |CurrentOutPtr|..|CurrentInPtr |CurrentOutPtr
//! ------------------------------------------------------------------------
//! OutPtr(buf2)|CurrentOutPtr|CurrentInPtr |..|CurrentOutPtr|CurrentInPtr
//! ------------------------------------------------------------------------
//! Result Buf  |    buf1     |    buf2       |..|        buf1    |    buf2
//! ------------------------------------------------------------------------
//!可见CurrentInPtr和CurrentOutPtr是交替的，这就是FFT中用两个缓冲区buf1和buf2进行乒乓操作，当阶数stage为奇数时结果存储在buf1中，当阶数stage为偶数时结果存储在buf2中，
//!  -# FFTSize must be a power of 2 (32, 64, 128, etc)
//!  -# FFTSize must be greater or equal to 32
//!  -# FFTStages must be log2(FFTSize)
//!  -# InPtr, OutPtr, CoefPtr, CurrentInPtr, CurrentOutPtr are FFTSize*2 in 
//!     length.
//!
int16_t main(void)
{
    // Locals
    uint16_t i;

#ifdef FLASH
    EALLOW;
    Flash0EccRegs.ECC_ENABLE.bit.ENABLE = 0;
    memcpy((uint32_t *)&RamfuncsRunStart, (uint32_t *)&RamfuncsLoadStart,
            (uint32_t)&RamfuncsLoadSize );
    FPU_initFlash();
#ifdef USE_TABLES//如果用到预先生成的旋转因子表，需要进行下列操作
    memcpy((uint32_t *)&FFTTwiddlesRunStart, (uint32_t *)&FFTTwiddlesLoadStart,
            (uint32_t)&FFTTwiddlesLoadSize );
#endif //USE_TABLES
#endif //FLASH
    
    FPU_initSystemClocks();
    
    FPU_initEpie();
    
    // Clear input buffers:
    for(i=0; i < (CFFT_SIZE*2); i=i+2){
        CFFTin1Buff[i] = 0.0f;
        CFFTin1Buff[i+1] = 0.0f;
        CFFTin2Buff[i] = 0.0f;
        CFFTin2Buff[i+1] = 0.0f;
        CFFToutBuff[i] = 0.0f;
        CFFToutBuff[i+1] = 0.0f;
    }
    
    // Generate sample waveforms://虚实交替存储
    //       CFFTin1Buff[0] = real[0]
    //       CFFTin1Buff[1] = imag[0]
    //       CFFTin1Buff[2] = real[1]
    //       厖�
    //       CFFTin1Buff[N] = real[N/2]
    //       CFFTin1Buff[N+1] = imag[N/2]
    //       厖�
    //       CFFTin1Buff[2N-3] = imag[N-2]
    //       CFFTin1Buff[2N-2] = real[N-1]
    //       CFFTin1Buff[2N-1] = imag[N-1]
    Rad = 0.0f;
    for(i=0; i < (CFFT_SIZE*2); i=i+2){
        CFFTin1Buff[i]   = sin(Rad) + cos(Rad*2.3567);       // Real Part
        CFFTin1Buff[i+1] = cos(Rad*8.345) + sin(Rad*5.789);  // Imaginary Part
        
        CFFTin2Buff[i]   = CFFTin1Buff[i];          // Not used in calculation
        CFFTin2Buff[i+1] = CFFTin1Buff[i+1];        // Not used in calculation
        Rad = Rad + RadStep;
    }

    //
    // Off-Place Algorithm异地算法
    //
    // Note: In this version, CFFTin1Buff and CFFToutBuff are used in
    //       ping-pong fashion. The input data is first stored in CFFTin1Buff
    //       where the FFT, including bit reversed ordering, is initially done.
    //       At each successive stage of the FFT the cfft.CurrentInPtr pointer 
    //       will point to the buffer that is the input for that stage. In this
    //       manner the "CurrentInPtr" and "CurrentOutPtr" are exchanged at the
    //       start of each stage. Depending on the number of FFT stages, the 
    //       final output will be in either CFFTin1Buff (#Stages is odd) or 
    //       CFFToutBuff (#stages is even).
    //注意：在此版本中，CFFTin1Buff和CFFToutBuff以乒乓球方式使用。输入数据首先存储在CFFTin1Buff中，其中FFT（包括位反向排序）最初完成。
//在FFT的每个连续阶段，CFFT。CurrentInPtr 指针将指向作为该阶段输入的缓冲区。通过这种方式，“CurrentInPtr”和“CurrentOutPtr”在每个阶段开始时交换。根据FFT级的数量，最终输出将以CFFTin1Buff（#Stages为奇数）或CFFToutBuff（#stages为偶数）为单位。
    //Input/output or middle stage of ping-pong buffer
    hnd_cfft->InPtr   = CFFTin1Buff;  
    //Output or middle stage of ping-pong buffer
    hnd_cfft->OutPtr  = CFFToutBuff;  
    hnd_cfft->Stages  = CFFT_STAGES;  // FFT stages
    hnd_cfft->FFTSize = CFFT_SIZE;    // FFT size
#ifdef USE_TABLES//用预先生成的旋转因子表CFFT_f32_twiddleFactors
    hnd_cfft->CoefPtr = CFFT_f32_twiddleFactors;  //Twiddle factor table
#else//不用预先生成的旋转因子表而是用 CFFT_f32_sincostable(hnd_cfft)自己生成旋转因子表
    hnd_cfft->CoefPtr = CFFTF32Coef;  //Twiddle factor table
    CFFT_f32_sincostable(hnd_cfft);   // Calculate twiddle factor
#endif //USE_TABLES
    
    //=========================================================================
    // CFFT result:
    //     CurrentInPtr[0] = real[0]
    //     CurrentInPtr[1] = imag[0]
    //     CurrentInPtr[2] = real[1]
    //     厖�
    //     CurrentInPtr[N] = real[N/2]
    //     CurrentInPtr[N+1] = imag[N/2]
    //     厖�
    //     CurrentInPtr[2N-3] = imag[N-2]
    //     CurrentInPtr[2N-2] = real[N-1]
    //     CurrentInPtr[2N-1] = imag[N-1]
    //
    //=========================================================================
#ifdef USE_TABLES
    CFFT_f32t(hnd_cfft);                   // Calculate FFT//计算FFT
#else
    CFFT_f32(hnd_cfft);                    // Calculate FFT
#endif //USE_TABLES
    
    // Check the output//CFFTgoldenOut是matlab计算结果的存储区、hnd_cfft->CurrentInPtr是DSP计算结果的存储区、
    for(i = 0; i < 2*CFFT_SIZE; i++){
        if(fabs(CFFTgoldenOut[i] - hnd_cfft->CurrentInPtr[i]) <= EPSILON){
            pass++;
        }else{
            fail++;
        }
    }
    //
    // Note: The input buffer for the magnitude calculation is pointed to by
    //       cfft.CurrentInPtr, while the output is stored in the memory 
    //       pointed to by cfft.CurrentOutPtr. If the user does not changed
    //       the value of cfft.CurrentOutPtr after calling the magnitude 
    //       calculation function, the output buffer is overwritten when the
    //       phase calculation function is called
    //幅度计算的输入缓冲区由 cfft.CurrentInPtr指向，而输出存储在cfft.CurrentOutPtr指向的内存中。如果用户未更改 cfft 的值。CurrentOutPtr 调用幅度计算函数后，调用相位计算函数时覆盖输出缓冲器
    //       If number of Stages is ODD, 如果FFT阶数为奇数
    //           currentInPtr=CFFTin1Buff, currentOutPtr=CFFToutBuff
    //       If number of Stages is EVEN, 
    //           currentInPtr=CFFToutBuff, currentOutPtr=CFFTin1Buff
    //

    // Calculate Magnitude:
#ifdef __TMS320C28XX_TMU__ //defined when --tmu_support=tmu0 in the project properties
                                         //project properties查看tmu_support
    // Calculate magnitude, result stored in CurrentOutPtr
    CFFT_f32_mag_TMU0(hnd_cfft);     
#else
    // Calculate magnitude, result stored in CurrentOutPtr    
    CFFT_f32_mag(hnd_cfft);          
#endif
    // Check the output
    for(i = 0; i < CFFT_SIZE; i++){
        if(fabs(CFFTgoldenMagnitude[i] - hnd_cfft->CurrentOutPtr[i]) <= 
                EPSILON){
            pass++;
        }else{
            fail++;
        }
    }
    // Magnitude Result:                       
    //       CurrentOutPtr[0] = Mag[0]           
    //       CurrentOutPtr[1] = Mag[1]           
    //       CurrentOutPtr[2] = Mag[2]           
    //       厖�                                 
    //       CurrentOutPtr[N-1] = Mag[N-1]       
                                           
    // Calculate Phase:
    // To avoid overwriting the magnitude, change the output buffer for 
    // the phase()
    hnd_cfft->CurrentOutPtr=CFFTin2Buff; 
                                         
#ifdef __TMS320C28XX_TMU__ //defined when --tmu_support=tmu0 in the project 
                           //properties
    // Calculate phase, result stored in CurrentOutPtr
    CFFT_f32_phase_TMU0(hnd_cfft);        
#else
    // Calculate phase, result stored in CurrentOutPtr    
    CFFT_f32_phase(hnd_cfft);              
#endif
    // Check the output
    for(i = 0; i < CFFT_SIZE; i++){
        if(fabs(CFFTgoldenPhase[i] - hnd_cfft->CurrentOutPtr[i]) <= EPSILON){
            pass++;
        }else{
            fail++;
        }
    }
    // Phase Result:
    //       CurrentOutPtr[0] = Phase[0]
    //       CurrentOutPtr[1] = Phase[1]
    //       CurrentOutPtr[2] = Phase[2]
    //       厖�
    //       CurrentOutPtr[N-1] = Phase[N-1]

    // End of test
    done();
    // Execution never reaches this point
    return 1;
}

以上是CFFT_f32_sincostable、CFFT_f32t(hnd_cfft)、CFFT_f32(hnd_cfft)、CFFT_f32_mag_TMU0(hnd_cfft)、CFFT_f32_mag(hnd_cfft)、CFFT_f32_phase_TMU0(hnd_cfft)、CFFT_f32_phase(hnd_cfft)几个库函数用法的例程
例程中运行的结果与例程的注释相符合

假设起初buf1为输入缓冲区，其地址由InPtr指向，buf2为输出缓冲区，其地址由OutPtr指向。则调用CFFT_f32t(hnd_cfft)或CFFT_f32(hnd_cfft)函数后：
如果FFT阶数为奇数，FFT的结果保存在输入缓冲区buf1中，其地址由InPtr和CurrentInPtr指向。
如果FFT阶数为偶数，FFT的结果保存在输入缓冲区buf2中，其地址由OutPtr和CurrentInPtr指向。

5.数据对齐

在cmd文件中实现了数据对齐
FFT阶数为7，所以做了N=2^7=128点的FFT，对于4N对齐应当使CFFT_ALIGNMENT=4*N=512。

6.关于#ifdef __ TMS320C28XX_TMU __

可以project properties查看tmu_support

在其中选择tmu0

结论：

假设起初buf1为输入缓冲区，其地址由InPtr指向，buf2为输出缓冲区，其地址由OutPtr指向。则调用CFFT_f32t(hnd_cfft)或CFFT_f32(hnd_cfft)函数后：
如果FFT阶数为奇数，FFT的结果保存在输入缓冲区buf1中，其地址由InPtr和CurrentInPtr指向。
如果FFT阶数为偶数，FFT的结果保存在输入缓冲区buf2中，其地址由OutPtr和CurrentInPtr指向。

---

于 2023-05-15 第一次整理编写

学习时整理，不当之处烦请指正
码字不易，留个赞再走吧