STM32单片机芯片与内部114 DSP-变换运算实数复数 FFT IFFT 不限制点数

cr4_fft_1024_stm32	函数用于实现1024点数据的FFT计算
cr4_fft_256_stm32	函数用于实现256点数据的FFT计算
cr4_fft_64_stm32	函数用于实现64点数据的FFT计算

2、计算幅频响应

uint32_t input[1024], output[1024], Mag[1024];/* 输入，输出和幅值 */
float32_t Phase[1024]; /* 相位*/
/*
*********************************************************************************************************
* 函 数 名: PowerMag
* 功能说明: 求模值
* 形 参： _usFFTPoints FFT点数
* 返 回 值: 无
*********************************************************************************************************
*/
void PowerMag(uint16_t _usFFTPoints)
{
    int16_t lX,lY;
    uint16_t i;
    float32_t mag;
    /* 计算幅值 */
    for (i=0; i < _usFFTPoints; i++)
    {
        lX= (output[i]<<16)>>16; /* 实部*/
        lY= (output[i]>> 16); /* 虚部 */
        arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
        Mag[i]= mag*2; /* 求模后乘以2才是实际模值，直流分量不需要乘2 */
    }
    /* 由于上面多乘了2，所以这里直流分量要除以2 */
    Mag[0]= Mag[0]>>1;
}

3、计算相频响应

/*
*********************************************************************************************************
* 函 数 名: Power_Phase_Radians
* 功能说明: 求相位
* 形 参： _usFFTPoints FFT点数， uiCmpValue 阀值
* 返 回 值: 无
*********************************************************************************************************
*/
void Power_Phase_Radians(uint16_t _usFFTPoints, uint32_t _uiCmpValue)
{
    int16_t lX, lY;
    uint16_t i;
    float32_t phase;
    float32_t mag;
    for (i=0; i <_usFFTPoints; i++)
    {
        lX= (output[i]<<16)>>16; /* 实部 */
        lY= (output[i] >> 16); /* 虚部 */
        phase = atan2(lY, lX); /* atan2求解的结果范围是(-pi, pi], 弧度制 */
        arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
        if(_uiCmpValue > mag)
        {
            Phase[i] = 0;
        }
        else
        {
            Phase[i] = phase* 180.0f/3.1415926f; /* 将求解的结果由弧度转换为角度 */
        }
}

二、复数浮点 FFT、IFFT（支持单精度和双精度）

新版 DSP 库浮点 FFT 推荐使用混合基函数 arm_cfft_f32，而基 2 函数 arm_cfft_radix2_f32 和基 4函数 arm_cfft_radix4_f32 将废弃。 ARM 说明如下：

DSP 库的早期发行版提供了单独的 radix-2 和 radix-4 对浮点数据进行运算的算法。这些功能仍然提供，但已弃用。相比新版函数，老版的功能较慢且通用性较低

1、基础支持

当前复数 FFT 函数支持三种数据类型，分别是浮点，定点 Q31 和 Q15。这些 FFT 函数有一个共同的特点，就是用于输入信号的缓冲，在转化结束后用来存储输出结果。这样做的好处是节省了 RAM 空间，不需要为输入和输出结果分别设置缓存。由于是复数 FFT，所以输入和输出缓存要存储实部和虚部。存储顺序如下： {real[0], imag[0], real[1], imag[1],………………} ，在使用中切记不要搞错。

2、单精度函数 arm_cfft_f32

支持16、32、64、128、256、512、1024、2048、4096点单精度复数FFT、IFFT。

3、双精度函数 arm_cfft_f64

支持16、32、64、128、256、512、1024、2048、4096点双精度复数FFT、IFFT

三、实数浮点 FFT（支持单精度和双精度）

CMSIS DSP 库里面包含一个专门用于计算实数序列的 FFT 库，很多情况下，用户只需要计算实数序列即可。计算同样点数 FFT 的实数序列要比计算同样点数的虚数序列有速度上的优势。

快速的 rfft 算法是基于混合基 cfft 算法实现的。

1、基础支持

一个 N 点的实数序列 FFT 正变换采用下面的步骤实现：

由上面的框图可以看出，实数序列的 FFT 是先计算 N/2 个实数的 CFFT，然后再重塑数据进行处理从而获得半个 FFT 频谱即可（利用了 FFT 变换后频谱的对称性）。

一个 N 点的实数序列 FFT 逆变换采用下面的步骤实现：

实数 FFT 支持浮点， Q31 和 Q15 三种数据类型。

2、单精度函数 arm_rfft_fast_f32

支持32、64、128、256、512、1024、2048、4096点单精度实数FFT、IFFT。

3、双精度函数 arm_rfft_fast_f64

支持32、64、128、256、512、1024、2048、4096点双精度实数FFT、IFFT

四、不限制点数 FFT 实现

可以看到前面的由于 ARM DSP 库限制最大只能 4096 点，有点不够用，所以整理了个更大点数的。不限制点数，满足 2^n 即可， n 最小值 4，即 16 个点的 FFT，而最大值不限。

此 FFT 算法没有再使用 ARM DSP 库，重新实现了一个。

1、函数 InitTableFFT

这个函数用于 FFT 计算过程中需要用到的正弦和余弦表。对于 8192 点和 16384 点已经专门制作了数值表，存到内部 Flash，其它点数继续使用的 RAM 空间。

/*
*********************************************************************************************************
* 函 数 名: 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

2、函数 cfft

这个函数用于复数 FFT 变换，如果需要实数则直接把虚部设置为0即可。

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) */
}