微波EDA网，见证研发工程师的成长！
首页 > 研发问答 > 嵌入式设计讨论 > MCU和单片机设计讨论 > 第32章 实数FFT的实现

第32章 实数FFT的实现

时间：10-02 整理：3721RD 点击：
第32章  实数FFT的实现

    本章主要讲解实数的浮点和定点Q31，Q15的实现。关于这部分的知识点和函数的计算结果上，官方的文档有一些小错误，在章节中会跟大家详细讲述，还有一个要注意的问题，调用实数FFT函数一定要使用CMSIS-DSP V1.4.4及其以上版本，以前的版本有bug。

    本章节使用的复数FFT函数来自ARM官方库的TransformFunctions部分

    32.1 复数FFT

    32.2 复数FFT-基2算法

    32.3 复数FFT-基4算法

    32.4 总结


32.1  实数FFT
32.1.1  描述

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

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

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



                              

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

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



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


32.2  实数FFT
32.2.1  arm_rfft_fast_f32

函数定义如下：

void arm_rfft_fast_f32(

arm_rfft_fast_instance_f32 * S,

float32_t * p, float32_t * pOut,

uint8_t ifftFlag)

参数定义：

[in] *S          points to anarm_rfft_fast_instance_f32 structure.

[in] *p          points to the inputbuffer.

[in] *pOut      points to the outputbuffer.

[in] ifftFlag     RFFT if flag is 0,RIFFT if flag is 1

注意事项：

结构arm_rfft_fast_instance_f32的定义如下（在文件arm_math.h文件）：

typedef struct

{

   arm_cfft_instance_f32 Sint;     /**< Internal CFFT structure. */

   uint16_t fftLenRFFT;            /**<length of the real sequence */

       float32_t* pTwiddleRFFT;          /**< Twiddle factors real stage  */

} arm_rfft_fast_instance_f32 ;

    下面通过在开发板上运行函数arm_rfft_fast_f32和arm_cfft_f32计算幅频响应，然后将相应的频率响应结果在Matlab上面绘制出来。

  1. /*
  2. *********************************************************************************************************
  3. *    函 数 名: arm_rfft_fast_f32_app
  4. *    功能说明: 调用函数arm_rfft_fast_f32计算1024点实数序列的幅频响应并跟使用函数arm_cfft_f32计算结果做对比
  5. *    形    参：无
  6. *    返 回 值: 无
  7. *********************************************************************************************************
  8. */
  9. static void arm_rfft_fast_f32_app(void)
  10. {
  11.      uint16_t i;
  12.      arm_rfft_fast_instance_f32 S;
  13.    
  14.      /* 实数序列FFT长度 */
  15.      fftSize = 1024;
  16.      /* 正变换 */
  17.     ifftFlag = 0;
  18.    
  19.      /* 初始化结构体S中的参数 */
  20.      arm_rfft_fast_init_f32(&S, fftSize);
  21.    
  22.      /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
  23.      for(i=0; i<1024; i++)
  24.      {
  25.          /* 50Hz正弦波，采样率1KHz */
  26.          testInput_f32_10khz[i] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;
  27.      }
  28.    
  29.      /* 1024点实序列快速FFT */
  30.      arm_rfft_fast_f32(&S, testInput_f32_10khz, testOutput_f32_10khz, ifftFlag);
  31.    
  32.      /* 为了方便跟函数arm_cfft_f32计算的结果做对比，这里求解了1024组模值，实际函数arm_rfft_fast_f32
  33.         只求解出了512组
  34.      */
  35.      arm_cmplx_mag_f32(testOutput_f32_10khz, testOutput, fftSize);

  37.      /* 串口打印求解的模值 */
  38.      for(i=0; i<fftSize; i++)
  39.      {
  40.          printf("%f\r\n", testOutput[i]);
  41.      }

  43.      printf("****************************分割线***************************************\r\n");

  45.      for(i=0; i<1024; i++)
  46.      {
  47.          /* 虚部全部置零 */
  48.          testInput_f32_10khz[i*2+1] = 0;
  49.         
  50.          /* 50Hz正弦波，采样率1KHz ,作为实部 */
  51.          testInput_f32_10khz[i*2] = 1.2f*arm_sin_f32(2*3.1415926f*50*i/1000)+1;
  52.      }
  53.    
  54.      arm_cfft_f32(&arm_cfft_sR_f32_len1024, testInput_f32_10khz, ifftFlag, doBitReverse);
  55.    
  56.      /* 求解模值  */
  57.      arm_cmplx_mag_f32(testInput_f32_10khz, testOutput, fftSize);

  59.      /* 串口打印求解的模值 */
  60.      for(i=0; i<fftSize; i++)
  61.      {
  62.          printf("%f\r\n", testOutput[i]);
  63.      }
  64.    
  65. }

复制代码

运行如上函数可以通过串口打印出函数arm_rfft_fast_f32和arm_cfft_f32计算的幅频模值，下面通过Matlab绘制波形来对比这两种模值。

    对比前需要先将串口打印出的两组数据加载到Matlab中，arm_rfft_fast_f32的计算结果起名signal，arm_cfft_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

  1. Fs = 1000;              % 采样率
  2. N  = 1024;            % 采样点数

  4. n  = 0:N-1;            % 采样序列
  5. f = n * Fs / N;          %真实的频率

  7. subplot(3,1,1);
  8. plot(f,  signal);        %绘制RFFT结果
  9. title('实数FFT');
  10. xlabel('时间');
  11. ylabel('幅值');

  13. subplot(3,1,2);
  14. plot(f,  sampledata);   %CFFT结果
  15. title('复数FFT');
  16. xlabel('时间');
  17. ylabel('幅值');

复制代码

Matlab运行结果如下：



                              

从上面的前512点对比中，我们可以看出两者的计算结果是相符的。这里有一点要特别注意，官方文档中对于函数arm_rfft_fast_f32输出结果的实部，虚部排列顺序说明是错误的。函数arm_rfft_fast_f32的输出结果仍然是实部，虚部，实部，虚部….. 依次排列下去。

    函数arm_rfft_fast_f32在计算直流分量（也就是频率为0的值）的虚部上是有错误的。关于这点大家可以将实际的实部和虚部输出结果打印出来做对比，但差别很小，基本可以忽略。


32.2.2  arm_rfft_q15

函数定义如下：

void arm_rfft_q15(

    const arm_rfft_instance_q15 * S,

    q15_t * pSrc,

    q15_t * pDst)

参数定义：

[in]  *S     points to an instance of the Q15 RFFT/RIFFTstructure.   

[in]  *pSrc  pointsto the input buffer.   

[out] *pDst  points to the output buffer.   

return       none.

注意事项：

结构arm_rfft_instance_q15的定义如下（在文件arm_math.h文件）：

typedef struct

{

uint32_t fftLenReal;   

uint8_t ifftFlagR;   

   uint8_t bitReverseFlagR;               

   uint32_t twidCoefRModifier;   

   q15_t *pTwiddleAReal;                  

   q15_t *pTwiddleBReal;                    

   const arm_cfft_instance_q15 *pCfft;      

} arm_rfft_instance_q15;

    下面通过在开发板上运行函数arm_rfft_q15和arm_cfft_f32计算幅频响应，然后将相应的频率响应结果在Matlab上面绘制出来。

  1. /*
  2. *********************************************************************************************************
  3. *    函 数 名: arm_rfft_q15_app
  4. *    功能说明: 调用函数arm_rfft_q15计算1024点实数序列的幅频响应并跟使用函数arm_cfft_f32计算的结果做对比。
  5. *    形    参：无
  6. *    返 回 值: 无
  7. *********************************************************************************************************
  8. */
  9. static void arm_rfft_q15_app(void)
  10. {
  11.      uint16_t i,j;
  12.      arm_rfft_instance_q15 S;
  13.    
  14.      /* 实数序列FFT长度 */
  15.      fftSize = 1024;
  16.      /* 正变换 */
  17.     ifftFlag = 0;
  18.      /* 码位倒序 */
  19.     doBitReverse = 1;
  20.    
  21.      /* 初始化结构体S */
  22.      arm_rfft_init_q15(&S, fftSize, ifftFlag, doBitReverse);
  23.    
  24.      /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
  25.      for(i=0; i<1024; i++)
  26.      {
  27.          /* 51.2Hz正弦波，采样率1024Hz。
  28.             arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
  29.             32768 / 20 = 1638.4
  30.          */
  31.          j = i % 20;
  32.          testInput_q15_50hz[i] = arm_sin_q15(1638*j);
  33.      }
  34.    
  35.      /* 1024点实序列快速FFT */
  36.      arm_rfft_q15(&S, testInput_q15_50hz, testOutput_q15_50hz);
  37.    
  38.      /* 由于输出结果的格式是Q5，所以这里将定点数转换为浮点数 */
  39.      for(i = 0; i < fftSize; i++)
  40.      {
  41.          testOutput_f32_10khz[i] = (float32_t)testOutput_q15_50hz[i]/32;
  42.      }
  43.    
  44.      /* 为了方便对比，这里求解了1024组复数，实际上面的变化只有512组
  45.         实际函数arm_rfft_q15只求解出了512组  */
  46.      arm_cmplx_mag_f32(testOutput_f32_10khz, testOutput, fftSize);

  48.      /* 串口打印求解的模值 */
  49.      for(i=0; i<fftSize; i++)
  50.      {
  51.          printf("%f\r\n", testOutput[i]);
  52.      }
  53.    
  54.      printf("****************************分割线***************************************\r\n");
  55.    
  56.      for(i=0; i<1024; i++)
  57.      {
  58.          /* 51.2Hz正弦波，采样率1024Hz。
  59.             arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
  60.             32768 / 20 = 1638.4
  61.          */
  62.          j = i % 20;
  63.          testInput_f32_10khz[i*2] = (float32_t) arm_sin_q15(1638*j)/32768;
  64.          /* 虚部全部置零 */
  65.          testInput_f32_10khz[i*2+1] = 0;
  66.         
  67.      }
  68.    
  69.      arm_cfft_f32(&arm_cfft_sR_f32_len1024, testInput_f32_10khz, ifftFlag, doBitReverse);
  70.    
  71.      /* 求解模值  */
  72.      arm_cmplx_mag_f32(testInput_f32_10khz, testOutput, fftSize);
  73.    
  74.      /* 串口打印求解的模值 */
  75.      for(i=0; i<fftSize; i++)
  76.      {
  77.          printf("%f\r\n", testOutput[i]);
  78.      }
  79. }

复制代码

运行如上函数可以通过串口打印出函数arm_rfft_q15和arm_cfft_f32计算的幅频模值，下面通过Matlab绘制波形来对比这两种模值。

     对比前需要先将串口打印出的两组数据加载到Matlab中，arm_rfft_q15的计算结果起名signal，arm_cfft_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

  1. Fs = 1000;             % 采样率
  2. N  = 1024;            % 采样点数

  4. n  = 0:N-1;            % 采样序列
  5. f = n * Fs / N;          %真实的频率

  7. subplot(3,1,1);
  8. plot(f,  signal);        %绘制RFFT结果
  9. title('实数FFT');
  10. xlabel('时间');
  11. ylabel('幅值');

  13. subplot(3,1,2);
  14. plot(f,  sampledata);   %CFFT结果
  15. title('复数FFT');
  16. xlabel('时间');
  17. ylabel('幅值');

复制代码

Matlab运行结果如下：



从上面的前512点对比中，我们可以看出两者的计算结果是相符的。这里有一点要特别注意，官方文档中对于函数arm_rfft_q31输出结果的实部，虚部排列顺序说明是错误的。函数arm_rfft_q31的输出结果仍然是实部，虚部，实部，虚部….. 依次排列下去。


32.2.3  arm_rfft_q31

函数定义如下：

void arm_rfft_q31(

    const arm_rfft_instance_q31 * S,

    q31_t * pSrc,

    q31_t * pDst)

参数定义：

[in]  *S     points to an instance of the Q31 RFFT/RIFFTstructure.   

[in]  *pSrc  pointsto the input buffer.   

[out] *pDst  points to the output buffer.   

return       none.

注意事项：

结构arm_rfft_instance_q31的定义如下（在文件arm_math.h文件）：

typedef struct

{

   uint32_t fftLenReal;                     

   uint8_t ifftFlagR;                          

   uint8_t bitReverseFlagR;                  

   uint32_t twidCoefRModifier;                 

   q31_t *pTwiddleAReal;                       

   q31_t *pTwiddleBReal;                       

   const arm_cfft_instance_q31 *pCfft;      

} arm_rfft_instance_q31;

    下面通过在开发板上运行函数arm_rfft_q31和arm_cfft_f32计算幅频响应，然后将相应的频率响应结果在Matlab上面绘制出来。

  1. /*
  2. *********************************************************************************************************
  3. *    函 数 名: arm_rfft_q31_app
  4. *    功能说明: 调用函数arm_rfft_q31计算1024点实数序列的幅频响应并跟使用函数arm_cfft_f32计算的结果做对比。
  5. *    形    参：无
  6. *    返 回 值: 无
  7. *********************************************************************************************************
  8. */
  9. static void arm_rfft_q31_app(void)
  10. {
  11.      uint16_t i,j;
  12.      arm_rfft_instance_q31 S;
  13.         
  14.      /* 实数序列FFT长度 */
  15.      fftSize = 1024;
  16.      /* 正变换 */
  17.     ifftFlag = 0;
  18.      /* 码位倒序 */
  19.     doBitReverse = 1;
  20.    
  21.      /* 初始化结构体S */
  22.      arm_rfft_init_q31(&S, fftSize, ifftFlag, doBitReverse);
  23.    
  24.      /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
  25.      for(i=0; i<1024; i++)
  26.      {
  27.         
  28.          /* 51.2Hz正弦波，采样率1024Hz。
  29.             arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
  30.             2^31 / 20 = 107374182.4
  31.          */
  32.          j = i % 20;
  33.          testInput_q31_50hz[i] = arm_sin_q31(107374182*j);
  34.      }
  35.    
  36.      /* 1024点实序列快速FFT */
  37.      arm_rfft_q31(&S, testInput_q31_50hz, testOutput_q31_50hz);
  38.    
  39.      /* 由于输出结果的格式是Q21，所以这里将定点数转换为浮点数 */
  40.      for(i = 0; i < fftSize; i++)
  41.      {
  42.          /* 输出的数据是11.21格式，2^21 = 4194304*/
  43.          testOutput_f32_10khz[i] = (float32_t)testOutput_q31_50hz[i]/2097152;
  44.      }
  45.    
  46.      /* 为了方便对比，这里求解了1024组复数，实际上面的变化只有512组
  47.         实际函数arm_rfft_q31只求解出了512组  */
  48.      arm_cmplx_mag_f32(testOutput_f32_10khz, testOutput, fftSize);
  49.    
  50.      /* 串口打印求解的模值 */
  51.      for(i=0; i<fftSize; i++)
  52.      {
  53.          printf("%f\r\n", testOutput[i]);
  54.      }
  55.    
  56.      printf("****************************分割线***************************************\r\n");
  57.    
  58.      for(i=0; i<1024; i++)
  59.      {
  60.          /* 51.2Hz正弦波，采样率1024Hz。
  61.             arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
  62.             2^31 / 20 = 107374182.4
  63.          */
  64.          j = i % 20;
  65.          testInput_f32_10khz[i*2] = (float32_t)arm_sin_q31(107374182*j)/2147483648;
  66.         
  67.          /* 虚部全部置零 */
  68.          testInput_f32_10khz[i*2+1] = 0;
  69.      }
  70.    
  71.      arm_cfft_f32(&arm_cfft_sR_f32_len1024, testInput_f32_10khz, ifftFlag, doBitReverse);
  72.    
  73.      /* 求解模值  */
  74.      arm_cmplx_mag_f32(testInput_f32_10khz, testOutput, fftSize);
  75.    
  76.      /* 串口打印求解的模值 */
  77.      for(i=0; i<fftSize; i++)
  78.      {
  79.          printf("%f\r\n", testOutput[i]);
  80.      }
  81. }

复制代码

运行如上函数可以通过串口打印出函数arm_rfft_q15和arm_cfft_f32计算的幅频模值，下面通过Matlab绘制波形来对比这两种模值。

     对比前需要先将串口打印出的两组数据加载到Matlab中，arm_rfft_q15的计算结果起名signal，arm_cfft_f32的计算结果起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

  1. Fs = 1000;             % 采样率
  2. N  = 1024;            % 采样点数

  4. n  = 0:N-1;            % 采样序列
  5. f = n * Fs / N;          %真实的频率

  7. subplot(3,1,1);
  8. plot(f,  signal);        %绘制RFFT结果
  9. title('实数FFT');
  10. xlabel('时间');
  11. ylabel('幅值');

  13. subplot(3,1,2);
  14. plot(f,  sampledata);   %CFFT结果
  15. title('复数FFT');
  16. xlabel('时间');
  17. ylabel('幅值');

复制代码

Matlab运行结果如下：



从上面的前512点对比中，我们可以看出两者的计算结果是相符的。这里有一点要特别注意，官方文档中对于函数arm_rfft_q31输出结果的实部，虚部排列顺序说明是错误的。函数arm_rfft_q31的输出结果仍然是实部，虚部，实部，虚部….. 依次排列下去。


32.3  总结

    使用实数FFT计算的时候要特别的注意本章节提到的几个错误点。


好资料。感谢分享

上一篇：急，NRF24L01+的使用
下一篇：五色音频谱源程序
栏目分类
热门文章

Copyright © 2017-2020 微波EDA网 版权所有

网站地图

Top