特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：链接
& R. e8 c4 B" s

第30章 复数FFT的实现

    本章主要讲解复数FFT的浮点和定点Q31，Q15的实现。

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

     30.1 复数FFT

    30.2 复数FFT-基2算法

    30.3 复数FFT-基4算法

    30.4 总结

1 [2 s; @" G* W+ q30.1 复数FFT30.1.1 描述

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

30.1.2 浮点

    浮点复数FFT使用了一个混合基数算法，通过多个基8与单个基2或基4算法实现。根据需要，该算法支持的长度[16，32，64，...，4096]和每个长度使用不同的旋转因子表。

    浮点复数FFT使用了标准的FFT定义，FFT正变换的输出结果会被放大fftLen倍数，计算FFT逆变换的时候会缩小到1/fftLen。这样就与教科书中的定义一致了。

    定义好的旋转因子和位反转表已经在头文件arm_const_structs.h中定义好了，调用浮点FFT函数arm_cfft_f32时，包含相应的头文件即可。比如：

       arm_cfft_f32(arm_cfft_sR_f32_len64, pSrc, 1, 1)

    上式就是计算一个64点的FFT逆变换包括位反转。数据结构arm_cfft_sR_f32_len64可以认为是常数，计算的过程中是不能修改的。同样是这种数据结构还能用于混合基的FFT正变换和逆变换。

    早期发布的浮点复数FFT函数版本包含基2和基4两种方法实现的，但是不推荐大家再使用了。现在全部用arm_cfft_f32代替了。

30.1.3 定点Q31和Q15

    定点库提供了基2和基4两种算法，基2算法支持的数据长度[16, 32, 64, ..., 4096]，基4算法支持的数据长度[16, 32, 64, ..., 4096]。一般情况下，建议使用基4算法，基4算法比基2算法执行速度要快一些。

    为了防止计算结果溢出，定点FFT每个蝶形运算的结果都要做放缩处理。对于基2算法，每次蝶形运算的结果要做0.5倍的放缩。对于基4算法，每次蝶形运算的结果要做0.25倍的放缩。FFT逆变换也要做相同的处理。相对于标准教科书式的FFT定义，定点FFT的计算结果放缩了1/fftLen（数据）倍。定点FFT的逆变也要做放缩处理，但是跟教科书式的FFT定义是相符的。

    每个FFT变换都需要一个单独的结构体，但结构体中的旋转因子和位反转表可以被重新使用。

    每个数据类型都有一个相关的初始化函数，初始化函数主要完成如下操作：

l 初始化结构体成员。

l 初始化旋转因子和位反转表指针。

    使用初始化函数是可选的。尽管如此，如果使用了初始化函数，那么结构体不能放在const data段，如果要放在const data段，应当按照如下方法进行初始化：

*arm_cfft_radix2_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};   

*arm_cfft_radix2_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};   

*arm_cfft_radix4_instance_q31 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};   

*arm_cfft_radix4_instance_q15 S = {fftLen, ifftFlag, bitReverseFlag, pTwiddle, pBitRevTable, twidCoefModifier, bitRevFactor};   

*arm_cfft_instance_f32 S = {fftLen, pTwiddle, pBitRevTable, bitRevLength};

    Q15和Q31 FFT使用了一个比较大的旋转因数和位反转表。这个表定义了一个最大长度的变换，表的子集可以用于短的变化。

30.1.4 arm_cfft_f32

函数定义如下：

    void arm_cfft_f32(

         const arm_cfft_instance_f32 * S,

         float32_t * p1,

         uint8_t ifftFlag,

         uint8_t bitReverseFlag)

参数定义：

     [in]   *S    points to an instance of the floating-point CFFT structure.  

     [in, out] *p1   points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.  

     [in] ifftFlag         flag that selects forward (ifftFlag=0) or inverse (ifftFlag=1) transform.  

     [in] bitReverseFlag  flag that enables (bitReverseFlag=1) or disables (bitReverseFlag=0) bit reversal of output.  

注意事项：

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

      typedef struct

      {

          uint16_t fftLen;                  

          const float32_t *pTwiddle;         

          const uint16_t *pBitRevTable;      

          uint16_t bitRevLength;            

      } arm_cfft_instance_f32;

0 s6 q8 ~# K2 a" [& H8 k6 L6 ^

    下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。

1. #define TEST_LENGTH_SAMPLES 2048
   ' J7 E' a( ]1 Y+ l* C* v3 D, Z  r
2. 
   
3. /* 输入和输出缓冲 */
   - p, ?. Q# D& W# x' H
4. static float32_t testOutput[TEST_LENGTH_SAMPLES/2];
   
5. static float32_t testInput_f32_10khz[TEST_LENGTH_SAMPLES];
   
6. 
   ( X* `6 M" x4 @1 h
7. /* 变量 */
   
8. uint32_t fftSize = 1024;
   1 L9 d) A" ?' L3 ?+ i
9. uint32_t ifftFlag = 0;
   * p$ J, v1 [" s" j) d
10. uint32_t doBitReverse = 1;
    $ i( |  Q# v* K: F4 u, [3 Q# `
11. 
    ; ~, u4 ~2 v  N% o; B# i
12. /*
    , D% A$ u7 r7 I$ c2 V4 B5 g; p
13. *********************************************************************************************************
    
14. *        函 数 名: arm_cfft_f32_app
    8 V  L. C' H" d
15. *        功能说明: 调用函数arm_cfft_f32_app计算幅频。
    - |6 Q' M/ B2 O
16. *        形    参：无
    
17. *        返 回 值: 无
    1 R" {) w2 p5 g; h
18. *********************************************************************************************************
    1 z6 \9 A2 ?4 x; k- [
19. */
    ; Z( K; U0 f5 w  H
20. void arm_cfft_f32_app(void)
    
21. {
    
22. uint16_t i;
    0 L( }  ~* K- e
23. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    , w7 Z$ V+ u2 q! m) K( O3 ~6 z% [
24. for(i=0; i<1024; i++)
    # a1 W& T' _) n3 R% p' A7 H, N
25. {
    9 A- Y/ Y6 F0 o5 T* ]. y
26. /* 虚部全部置零 */
    8 N5 P6 Y) g  i1 F
27. testInput_f32_10khz[i*2+1] = 0;
    " T& p$ G( E: v! P% u
28. /* 50Hz正弦波，采样率1KHz ,作为实部 */
    
29. testInput_f32_10khz[i*2] = arm_sin_f32(2*3.1415926f*50*i/1000);
    
30. }
    # o# O4 c' [) A6 c: A
31. /* CFFT变换 */
    ) V: |: A. f* l
32. arm_cfft_f32(&arm_cfft_sR_f32_len1024, testInput_f32_10khz, ifftFlag, doBitReverse);
    
33. 
    
34. /* 求解模值  */
    
35. arm_cmplx_mag_f32(testInput_f32_10khz, testOutput, fftSize);
    ! b+ U3 f% s0 G9 R
36. /* 串口打印求解的模值 */
    1 |. k4 W2 I& p( ]7 T
37. for(i=0; i<1024; i++)
    
38. {
    
39. printf("%f\r\n", testOutput[i]);
    * ~, ^/ A% Q; Z; a
40. }
    * Z: H' x" \# R8 o2 r3 @' l  m7 T
41. 
    $ y9 N+ B- b* D2 }
42. }

复制代码

运行如上函数可以通过串口打印出计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_f32计算的模值做对比。

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

Fs = 1000;                % 采样率

N  = 1024;               % 采样点数

n  = 0:N-1;               % 采样序列

t  = 0:1/Fs:1-1/Fs;        % 时间序列

f = n * Fs / N;             %真实的频率

0 ?, C: Y0 A3 O& f

x = sin(2*pi*50*t) ;  %原始信号

y = fft(x, N);        %对原始信号做FFT变换

: J' O, G2 Z; J& R- \

subplot(2,1,2);

plot(f, abs(y));       %绘制幅频相应曲线

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

4 R1 k$ k% Q2 _' c
- X+ ^1 j% n7 c  l8 {1 j, Z

subplot(2,1,1);

plot(f,  sampledata);       %绘制幅频相应曲线

title('复数FFT计算结果');

xlabel('频率');

ylabel('幅度');

8 M2 A3 W  I/ {, Z9 W; m; @8 O2 h( S
8 E2 I6 J& H0 w& ^

Matlab运行结果如下：

从上面的对比结果中可以看出，Matlab和函数arm_cfft_f32计算的结果基本是一直的。

9 N- O7 C. F" e" u3 t

30.2 复数FFT—基2算法
30.2.1 arm_cfft_radix2_f32

    此函数已经不推荐使用，后面的版本会被删除，故不做介绍。

30.2.2 arm_cfft_radix2_q31

函数定义如下：

    void arm_cfft_radix2_q31(

    const arm_cfft_radix2_instance_q31 * S,

    q31_t * pSrc)

参数定义：

    [in]      *S    points to an instance of the fixed-point CFFT/CIFFT structure.  

    [in, out] *pSrc  points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.

注意事项：

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

      typedef struct

      {

          uint16_t fftLen;               

          uint8_t ifftFlag;               

          uint8_t bitReverseFlag;         

          q31_t *pTwiddle;                 

          uint16_t *pBitRevTable;         

          uint16_t twidCoefModifier;      

          uint16_t bitRevFactor;         

      } arm_cfft_radix2_instance_q31;

    下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。

1. q31_t testInput_radix2_q31_50hz[TEST_LENGTH_SAMPLES];
   
2. q31_t testOutputQ31[TEST_LENGTH_SAMPLES/2];
   
3. 
   ( h) m8 m& A5 Y! ~2 m" J/ i8 N0 Z
4. /*
   
5. *********************************************************************************************************
   . i0 \) c' n- Z5 ^, K) E
6. *        函 数 名: arm_cfft_radix2_q31_app
   
7. *        功能说明: 调用函数arm_cfft_radix2_q31计算幅频。
   " q; ~* w  F$ e0 U  x& G
8. *        形    参：无
   3 |6 X1 z0 k) O' ^# ?" w2 R
9. *        返 回 值: 无
   
10. *********************************************************************************************************
    2 l. s/ p3 D3 y: W. i3 Y4 a
11. */
    1 m, L4 _( R, q1 U  p  O7 \' [# S
12. void arm_cfft_radix2_q31_app(void)
    
13. {
    
14. uint16_t i,j;
    
15. arm_cfft_radix2_instance_q31 S;
    2 A8 ]( b5 C# Y1 F1 q& e, A
16. fftSize = 1024;
    - R9 \- y( R9 ]& }% U
17.     ifftFlag = 0;
    , ?" T+ |/ w* t
18.     doBitReverse = 1;
    
19. /* 初始化结构S */
    
20. arm_cfft_radix2_init_q31(&S, fftSize, ifftFlag, doBitReverse);
    ; f) q# d, c4 z# @
21. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    
22. for(i=0; i<1024; i++)
    
23. {
    , }0 E7 E& S* H! x5 D1 C, q
24. testInput_radix2_q31_50hz[i*2+1] = 0;
    
25. /* 51.2Hz正弦波，采样率1024Hz。
    
26.    arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
    
27.    2^31 / 20 = 107374182.4
    
28. */
    : V4 G& f$ C* t/ `# P9 t
29. j = i % 20;
    2 _2 K8 ^( Z4 O8 i0 ], K! n
30.           testInput_radix2_q31_50hz[i*2] = arm_sin_q31(107374182*j);
    9 t7 V' F) T0 V
31.           printf("%d\r\n", testInput_radix2_q31_50hz[i*2]);
    , J: k$ A$ [6 I" g- i2 [  U. R# T
32. }
    
33. /* 输出结果分割线 */
    ! Y$ u& O6 p( V! f8 h' r# l1 ]( J
34. printf("**************************************************\r\n");
    
35. printf("**************************************************\r\n");
      ^- _  O8 I! O9 S2 u4 W+ f
36. /* 计算CFFT */
    
37. arm_cfft_radix2_q31(&S, testInput_radix2_q31_50hz);
    ' w8 _& s, U( D  f) u
38. /* 计算模值 */
    
39. arm_cmplx_mag_q31(testInput_radix2_q31_50hz, testOutputQ31, fftSize);
    # j4 C3 x* o3 }; J
40. /* 串口打印求解的模值 */
    % @, W, u" L; q6 j* p
41. for(i=0; i<1024; i++)
    
42. {
    , z/ d, e' d, G9 j/ [/ w9 J
43. printf("%d\r\n", testOutputQ31[i]);
    2 O: \3 y- h/ T4 N. _
44. }
    
45. }

复制代码

运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix2_q31计算的模值做对比。

    对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q31计算的模值起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

Fs = 1024;              % 采样率

N  = 1024;             % 采样点数

n  = 0:N-1;             % 采样序列

f = n * Fs / N;           %真实的频率

5 u' c: c2 H* |  ]
: F( A; O  Z: I; w3 @' }  f0 K2 y  w

y = fft(signal, N);        %对原始信号做FFT变换

% T$ G8 b7 C+ O6 t1 t. K
) S  p/ p4 Q  |1 N/ u* g

subplot(2,1,2);

plot(f, abs(y));           %绘制幅频相应曲线

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

subplot(2,1,1);

plot(f,  sampledata);      %绘制幅频相应曲线

title('复数FFT计算结果');

xlabel('频率');

ylabel('幅度');

Matlab运行结果如下：

从上面的对比结果中可以看出，Matlab和函数arm_cfft_radix2_q31计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q31对数据结果做了移位处理。

30.2.3 arm_cfft_radix2_q15

函数定义：

    void arm_cfft_radix2_q15(

        const  arm_cfft_radix2_instance_q15 * S,

        q15_t  * pSrc)

参数定义：

    [in]      *S    points to an instance of the fixed-point CFFT/CIFFT structure.  

    [in, out] *pSrc  points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.  

注意事项：

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

      typedef struct

      {

         uint16_t fftLen;               

          uint8_t ifftFlag;               

          uint8_t bitReverseFlag;         

          q15_t *pTwiddle;                    

      uint16_t *pBitRevTable;         

      uint16_t twidCoefModifier;      

          uint16_t bitRevFactor;         

      } arm_cfft_radix2_instance_q15;

8 c6 a$ G8 l6 d. l  d, `5 G/ [

    下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。

1. q15_t testInput_radix2_q15_50hz[TEST_LENGTH_SAMPLES];
   2 w7 Y' |' }5 S8 V2 t
2. q15_t testOutputQ15[TEST_LENGTH_SAMPLES/2];
   0 o$ Y$ m+ {: @. A2 m% r0 \
3. 
   
4. /*
   - i( [# w' J* ^3 k
5. *********************************************************************************************************
   
6. *        函 数 名: arm_cfft_radix2_q15_app
   
7. *        功能说明: 调用函数arm_cfft_radix2_q15计算幅频。
   ) C6 a; M& g# o7 }3 R: p( |
8. *        形    参：无
   8 L0 D1 D  b1 D! Q; M' S
9. *        返 回 值: 无
   
10. *********************************************************************************************************
    + v( B0 I0 m; {. X- j* V6 o
11. */
    1 ?3 f: s. ]& q5 v4 q9 D  }
12. void arm_cfft_radix2_q15_app(void)
    
13. {
    
14. uint16_t i,j;
    
15. arm_cfft_radix2_instance_q15 S;
    5 ?( j% @; P$ B- V! p
16. fftSize = 1024;
    ' D- y# Y0 W# w7 t
17.     ifftFlag = 0;
    
18.     doBitReverse = 1;
    9 G8 d; H9 L  ]( C! _7 t
19. /* 初始化结构S */
    
20. arm_cfft_radix2_init_q15(&S, fftSize, ifftFlag, doBitReverse);
    % Y& D# c7 i1 }8 V3 D2 v# @8 \
21. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    ) A) X& y9 R6 r4 l! y' [0 t
22. for(i=0; i<1024; i++)
    1 @' F0 C% j: s4 p, x, N
23. {
    5 r: Q6 x8 a, B" o
24. /* 虚部全部置0 */
    
25. testInput_radix2_q15_50hz[i*2+1] = 0;
    
26. /* 51.2Hz正弦波，采样率1024Hz。
    ' T# b6 R/ M  d3 q- }2 ?
27.    arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
    
28.    32768 / 20 = 1638.4
    2 p+ |- r9 q# s6 q9 h
29. */
    * f& G! V$ [, l
30. j = i % 20;
    
31.           testInput_radix2_q15_50hz[i*2] = arm_sin_q15(1638*j);
    0 {5 a$ i  z. u: _6 X" r4 N
32.           printf("%d\r\n", testInput_radix2_q15_50hz[i*2]);
    : ?& N# s0 E- _5 i
33. }
    $ N, W7 U1 g5 J9 ?( @$ h' x+ t
34. /* 输出结果分割线 */
    
35. printf("**************************************************\r\n");
    1 E9 E& n6 [) J
36. printf("**************************************************\r\n");
    - E: n. |& I% k  R( E0 q
37. /* 计算CFFT */
    
38. arm_cfft_radix2_q15(&S, testInput_radix2_q15_50hz);
    ! l1 r9 V/ J2 C& ~0 D& t
39. /* 计算模值 */
      o1 l# K) N! [) u) a
40. arm_cmplx_mag_q15(testInput_radix2_q15_50hz, testOutputQ15, fftSize);
    ' ~+ D& m+ l' ?
41. /* 串口打印求解的模值 */
    
42. for(i=0; i<1024; i++)
    5 c9 V( O$ u1 T& B! }6 |5 y
43. {
    4 p$ J- y9 {! r' e6 ?
44. printf("%d\r\n", testOutputQ15[i]);
    
45. }
    
46. }

复制代码

运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix2_q15计算的模值做对比。

    对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q15计算的模值起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

Fs = 1024;                     % 采样率

N  = 1024;                     % 采样点数

n  = 0:N-1;                     % 采样序列

f = n * Fs / N;                %真实的频率

y = fft(signal, N);           %对原始信号做FFT变换

2 G9 i! Q* ^. v* t; U+ m4 C5 n- ?

subplot(2,1,2);

plot(f, abs(y));           %绘制幅频相应曲线

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

' N9 c% J4 u; F9 [5 ?

subplot(2,1,1);

plot(f,  sampledata);      %绘制幅频相应曲线

title('复数FFT计算结果');

xlabel('频率');

ylabel('幅度');

9 G5 w5 z. e* Q: W! r$ H& x' M

Matlab运行结果如下：

从上面的对比结果中可以看出，Matlab和函数arm_cfft_radix2_q15计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q15对数据结果做了移位处理。

30.3 复数FFT—基4算法

    如果希望直接调用FFT程序计算IFFT，可以用下面的方法：

30.3.1 arm_cfft_radix4_f32

    此函数已经不推荐使用，后面的版本会被删除，故不做介绍。

30.3.2 arm_cfft_radix4_q31

函数定义如下：

    void arm_cfft_radix4_q31(

    const arm_cfft_radix4_instance_q31 * S,

    q31_t * pSrc)

参数定义：

    [in]      *S    points to an instance of the fixed-point CFFT/CIFFT structure.  

    [in, out] *pSrc  points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.

注意事项：

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

      typedef struct

      {

          uint16_t fftLen;               

          uint8_t ifftFlag;               

          uint8_t bitReverseFlag;         

          q31_t *pTwiddle;                 

          uint16_t *pBitRevTable;         

          uint16_t twidCoefModifier;      

          uint16_t bitRevFactor;         

       } arm_cfft_radix4_instance_q31;

2. 为了防止数据饱和，每次蝶形运行的结果都要除以2，故不同的长度的FFT运算的最终结果输出格式不同。具体信息如下：

    Q31 CFFT

    Q31 CIFFT

    下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。

1. q31_t testInput_radix4_q31_50hz[TEST_LENGTH_SAMPLES];
   . d- N* ^; `6 l! ?5 Y) q$ k
2. /*
   
3. *********************************************************************************************************
   
4. *        函 数 名: arm_cfft_radix4_q31_app
   
5. *        功能说明: 调用函数arm_cfft_radix4_q31_app计算幅频。
   5 y2 `) Z% Z' N( }( w' `7 a; h! Z
6. *        形    参：无
   
7. *        返 回 值: 无
   
8. *********************************************************************************************************
   
9. */
   ; A( V3 W9 [! c9 _1 ~3 T6 ~: B
10. void arm_cfft_radix4_q31_app(void)
    7 D. h& p, ?+ r9 ?
11. {
      b4 j) l% S4 ~& J/ x: y
12. uint16_t i,j;
    
13. arm_cfft_radix4_instance_q31 S;
    1 v6 y, @& _; e7 J, f
14. fftSize = 1024;
    
15.     ifftFlag = 0;
    
16.     doBitReverse = 1;
    
17. /* 初始化结构S */
    
18. arm_cfft_radix4_init_q31(&S, fftSize, ifftFlag, doBitReverse);
    ' J/ a' e3 N7 b1 G! B1 ^
19. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    4 ^! E/ ]: O) P4 p. @
20. for(i=0; i<1024; i++)
    
21. {
    : y2 S5 z# A/ w/ D1 p
22. testInput_radix4_q31_50hz[i*2+1] = 0;
    
23. /* 51.2Hz正弦波，采样率1024Hz。
    ! R4 F  ]% H- v; ~
24.    arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
    
25.    2^31 / 20 = 107374182.4
    
26. */
    ' u, v, z/ N$ o. N( v
27. j = i % 20;
    
28.           testInput_radix4_q31_50hz[i*2] = arm_sin_q31(107374182*j);
    : {; l& p9 L" X9 L; x- b; ^
29.           printf("%d\r\n", testInput_radix4_q31_50hz[i*2]);
    5 g/ Z/ Q! a% L3 C! }
30. }
    # u! B! ~9 u; J0 a- `5 n
31. /* 输出结果分割线 */
    3 |! s" h8 U6 r! x
32. printf("**************************************************\r\n");
    
33. printf("**************************************************\r\n");
    
34. /* 计算CFFT */
    + t$ c/ N( Y) y
35. arm_cfft_radix4_q31(&S, testInput_radix4_q31_50hz);
    
36. /* 计算模值 */
    
37. arm_cmplx_mag_q31(testInput_radix4_q31_50hz, testOutputQ31, fftSize);
    7 l: r8 @8 ~' y( B
38. /* 串口打印求解的模值 */
    
39. for(i=0; i<1024; i++)
    & f8 }% V' J$ ~' G
40. {
    ; I. H8 p8 S5 [' l
41. printf("%d\r\n", testOutputQ31[i]);
    2 O( _1 u. Z% e
42. }
    1 [1 ^4 J& L1 W% k
43. }

复制代码

运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix4_q31计算的模值做对比。

    对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q31计算的模值起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

Fs = 1024;              % 采样率

N  = 1024;             % 采样点数

n  = 0:N-1;             % 采样序列

f = n * Fs / N;           %真实的频率

9 L9 V/ N. b. T- c8 H

y = fft(signal, N);        %对原始信号做FFT变换

. a" h, ~5 q1 c8 r

subplot(2,1,2);

plot(f, abs(y));           %绘制幅频相应曲线

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

3 ?( `. M$ B& s

subplot(2,1,1);

plot(f,  sampledata);      %绘制幅频相应曲线

title('复数FFT计算结果');

xlabel('频率');

ylabel('幅度');

" `" w7 w. V6 y& z2 i0 L# x2 C

Matlab运行结果如下：

从上面的对比结果中可以看出，Matlb和函数arm_cfft_radix4_q31计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q31和arm_cfft_radix4_q31对数据结果做了移位处理。

30.3.3 arm_cfft_radix4_q15

函数定义：

    void arm_cfft_radix4_q15(

        const  arm_cfft_radix4_instance_q15 * S,

        q15_t  * pSrc)

参数定义：

    [in]      *S    points to an instance of the fixed-point CFFT/CIFFT structure.  

    [in, out] *pSrc  points to the complex data buffer of size <code>2*fftLen</code>. Processing occurs in-place.  

注意事项：

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

      typedef struct

      {

          uint16_t fftLen;               

          uint8_t ifftFlag;               

          uint8_t bitReverseFlag;         

          q15_t *pTwiddle;                    

      uint16_t *pBitRevTable;         

      uint16_t twidCoefModifier;      

          uint16_t bitRevFactor;         

      } arm_cfft_radix4_instance_q15;

2. 为了防止数据饱和，每次蝶形运行的结果都要除以2，故不同的长度的FFT运算的最终结果输出格式不同。具体信息如下：

    Q15 CFFT

    Q15 CIFFT

    下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。

1. q15_t testInput_radix4_q15_50hz[TEST_LENGTH_SAMPLES];
   3 t/ o& J' [. r9 y
2. /*
   
3. *********************************************************************************************************
   3 w4 n( n. H% x6 T/ c2 y  v" g8 G
4. *        函 数 名: arm_cfft_radix4_q15_app
   * G7 k  l* V& k5 ?; ]
5. *        功能说明: 调用函数arm_cfft_radix4_q15计算幅频。
   2 c: k$ R9 Q( |4 r" M
6. *        形    参：无
   9 P# Z# V2 m- x6 }" w! O- X
7. *        返 回 值: 无
   
8. *********************************************************************************************************
   9 T  N' W% |! [: f! }  i
9. */
   
10. void arm_cfft_radix4_q15_app(void)
    
11. {
    6 x' o, W1 Y7 L
12. uint16_t i,j;
    * O/ f1 o0 m5 ?% s
13. arm_cfft_radix4_instance_q15 S;
    
14. fftSize = 1024;
    
15.     ifftFlag = 0;
    
16.     doBitReverse = 1;
    ' c) L& n) F: I2 Q$ f6 J' D: W7 W
17. /* 初始化结构S */
    
18. arm_cfft_radix4_init_q15(&S, fftSize, ifftFlag, doBitReverse);
    
19. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    
20. for(i=0; i<1024; i++)
    
21. {
    
22. /* 虚部全部置0 */
    
23. testInput_radix4_q15_50hz[i*2+1] = 0;
    4 k* _4 C0 ]- D
24. /* 51.2Hz正弦波，采样率1024Hz。
    ' _6 y5 ?" i- B
25.    arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
    , v# G( g6 J+ x
26.    32768 / 20 = 1638.4
    
27. */
    
28. j = i % 20;
    ! N  u, ^, x3 s
29.           testInput_radix4_q15_50hz[i*2] = arm_sin_q15(1638*j);
    , B# L; b- D1 r8 X$ {
30.           printf("%d\r\n", testInput_radix4_q15_50hz[i*2]);
    
31. }
    $ e! e& Z' O) N( Q' o. ^' J
32. /* 输出结果分割线 */
    
33. printf("**************************************************\r\n");
    
34. printf("**************************************************\r\n");
    
35. /* 计算CFFT */
    2 m0 r7 a1 U) ^
36. arm_cfft_radix4_q15(&S, testInput_radix4_q15_50hz);
    
37. /* 计算模值 */
    
38. arm_cmplx_mag_q15(testInput_radix4_q15_50hz, testOutputQ15, fftSize);
    
39. /* 串口打印求解的模值 */
      J' X) _1 [: y9 l
40. for(i=0; i<1024; i++)
    5 U; @$ i$ K$ ?8 B' a7 N5 j$ }2 L
41. {
    2 \( g3 [, z' @9 W) e/ b
42. printf("%d\r\n", testOutputQ15[i]);
    - q, X. F9 v* S$ k% M5 ~" n' O
43. }
    
44. }

复制代码

运行如上函数可以通过串口打印出原始信号和计算的模值，下面我们就通过Matlab计算的模值跟arm_cfft_radix4_q15计算的模值做对比。

    对比前需要先将串口打印出的数据加载到Matlab中，原始信号起名signal，函数arm_cmplx_mag_q15计算的模值起名sampledata，加载方法在前面的教程中已经讲解过，这里不做赘述了。Matlab中运行的代码如下：

Fs = 1024;                     % 采样率

N  = 1024;                     % 采样点数

n  = 0:N-1;                     % 采样序列

f = n * Fs / N;                %真实的频率

& D0 E; e1 a' [; @- f

y = fft(signal, N);           %对原始信号做FFT变换

" F9 c3 {; C  {4 S
* P+ R5 q9 T4 Y  o' k

subplot(2,1,2);

plot(f, abs(y));           %绘制幅频相应曲线

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

subplot(2,1,1);

plot(f,  sampledata);      %绘制幅频相应曲线

title('复数FFT计算结果');

xlabel('频率');

ylabel('幅度');

, Q$ b& s" X) N4 B: X4 `+ @9 N

Matlab运行结果如下：

从上面的对比结果中可以看出，Matlab和函数arm_cfft_radix4_q15计算的频率点基本是一致的，而幅值大小不一样是因为调用函数arm_cmplx_mag_q15和arm_cfft_radix4_q15对数据结果做了移位处理。

1 }& u' V) W. l( y  T- F7 u

30.4 总结

    本章节内容比较多，有兴趣的可以深入了解源码的实现。

标记下学习

mark一下！