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

第30章 复数FFT的实现

' M* A, r9 O4 t6 Z

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

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

     30.1 复数FFT

    30.2 复数FFT-基2算法

    30.3 复数FFT-基4算法

    30.4 总结

2 }" y+ h* i( v3 R30.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;

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

1. #define TEST_LENGTH_SAMPLES 2048
   - C/ b: x  k4 t6 T! B, `
2. 
   
3. /* 输入和输出缓冲 */
   
4. static float32_t testOutput[TEST_LENGTH_SAMPLES/2];
   3 ]# h2 C3 _- x9 K! e- c  p
5. static float32_t testInput_f32_10khz[TEST_LENGTH_SAMPLES];
   
6. 
   
7. /* 变量 */
   
8. uint32_t fftSize = 1024;
   
9. uint32_t ifftFlag = 0;
   
10. uint32_t doBitReverse = 1;
    , v4 |5 y5 Y, Z8 @. l8 l
11. 
    6 u, d. `5 O) o5 j6 k7 N
12. /*
    
13. *********************************************************************************************************
    8 {" e  y+ |. l& H
14. *        函 数 名: arm_cfft_f32_app
    ; X6 t6 c+ y- B) g8 u5 R( }3 `& S
15. *        功能说明: 调用函数arm_cfft_f32_app计算幅频。
    1 L1 C" _2 G4 `) Y) [" w3 y
16. *        形    参：无
    + G+ e8 C3 m6 v) R* b& K+ o
17. *        返 回 值: 无
    6 j8 P) T: s& K. b8 `
18. *********************************************************************************************************
    
19. */
    
20. void arm_cfft_f32_app(void)
    & }, h4 M1 H' Q0 `2 ~  p/ M$ z' E$ }
21. {
    
22. uint16_t i;
    
23. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    
24. for(i=0; i<1024; i++)
    * s- X7 t/ H7 |
25. {
    5 i$ v5 r" m+ G/ i8 O; F9 h8 L
26. /* 虚部全部置零 */
    6 _1 W* R# `2 B) r. X/ t+ h
27. testInput_f32_10khz[i*2+1] = 0;
    
28. /* 50Hz正弦波，采样率1KHz ,作为实部 */
    7 R7 B/ m: w$ v* _
29. testInput_f32_10khz[i*2] = arm_sin_f32(2*3.1415926f*50*i/1000);
    
30. }
    
31. /* CFFT变换 */
    3 R6 K" b, V! f( M/ Z/ |% K. K
32. arm_cfft_f32(&arm_cfft_sR_f32_len1024, testInput_f32_10khz, ifftFlag, doBitReverse);
    1 y$ G* I0 \% S" `& ~6 S
33. 
    
34. /* 求解模值  */
    1 s- s  G6 E& T
35. arm_cmplx_mag_f32(testInput_f32_10khz, testOutput, fftSize);
    6 n& b; f5 l1 e2 M/ F7 j
36. /* 串口打印求解的模值 */
    : m( z/ u- W: E+ l2 j  h8 ?
37. for(i=0; i<1024; i++)
    
38. {
    
39. printf("%f\r\n", testOutput[i]);
    6 j/ |( w  c2 f- \# y- n8 E
40. }
    / i7 C# [, y) @; @- q, F
41. 
      P' n$ V; C; Q! [
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;             %真实的频率

& l" [# e# u5 e5 [/ u9 h# R

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

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

2 W3 y2 }; }3 s* |# n

subplot(2,1,2);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

) Z4 Y: ~9 Y$ y+ Q

subplot(2,1,1);

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

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

xlabel('频率');

ylabel('幅度');

& ?- p; C9 s9 t* `: l& |- Z0 a1 h

Matlab运行结果如下：

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

3 f8 ^  o1 J: I! m' S1 p

30.2 复数FFT—基2算法
+ i# h: F% D6 P6 A' C/ @30.2.1 arm_cfft_radix2_f32

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

+ g) a; S7 N* D30.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];
   1 y/ ~0 e/ Q/ f$ ]- o
2. q31_t testOutputQ31[TEST_LENGTH_SAMPLES/2];
   
3. 
   
4. /*
   
5. *********************************************************************************************************
   % F0 n% z! z7 A  w0 T& U) T4 U+ A/ k( R+ W
6. *        函 数 名: arm_cfft_radix2_q31_app
   
7. *        功能说明: 调用函数arm_cfft_radix2_q31计算幅频。
   
8. *        形    参：无
   
9. *        返 回 值: 无
   8 n! P* d) Y# u0 V" o& w% I
10. *********************************************************************************************************
    
11. */
    
12. void arm_cfft_radix2_q31_app(void)
    
13. {
    ! l8 K' |2 C$ u$ r# m% [- i% c. f, R
14. uint16_t i,j;
    ' P* q+ p2 a3 u( o3 g% v- i
15. arm_cfft_radix2_instance_q31 S;
    
16. fftSize = 1024;
    
17.     ifftFlag = 0;
    
18.     doBitReverse = 1;
    ' A3 v& `/ m5 F' X
19. /* 初始化结构S */
    % B# {% W8 O9 h7 E
20. arm_cfft_radix2_init_q31(&S, fftSize, ifftFlag, doBitReverse);
    7 ~8 b. h% L2 a5 ^
21. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    
22. for(i=0; i<1024; i++)
    
23. {
    7 Z5 _' O5 {0 y- R7 M# Z8 T8 o+ B
24. testInput_radix2_q31_50hz[i*2+1] = 0;
    
25. /* 51.2Hz正弦波，采样率1024Hz。
    " V! u* v3 s  s5 J" {
26.    arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
    
27.    2^31 / 20 = 107374182.4
    % h; q9 F3 g5 `7 B: b
28. */
    
29. j = i % 20;
    
30.           testInput_radix2_q31_50hz[i*2] = arm_sin_q31(107374182*j);
    
31.           printf("%d\r\n", testInput_radix2_q31_50hz[i*2]);
    
32. }
    , y* K1 W0 @8 o5 _+ l
33. /* 输出结果分割线 */
    ! i9 f+ h3 C, A, t- D" s1 X
34. printf("**************************************************\r\n");
    
35. printf("**************************************************\r\n");
    
36. /* 计算CFFT */
    8 ~3 h! `% P* b
37. arm_cfft_radix2_q31(&S, testInput_radix2_q31_50hz);
    
38. /* 计算模值 */
    
39. arm_cmplx_mag_q31(testInput_radix2_q31_50hz, testOutputQ31, fftSize);
    ' X+ ~: y% v7 k1 E8 |! Z0 }& t0 g/ g
40. /* 串口打印求解的模值 */
    
41. for(i=0; i<1024; i++)
    $ s+ V' K8 @0 i6 Q* m, x
42. {
    / B! a8 H9 J' n% Z" o
43. printf("%d\r\n", testOutputQ31[i]);
    * F. r% v. G6 o) D4 O
44. }
    8 N/ g( z6 Z5 A1 [- T# ~1 _
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;           %真实的频率

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

. W4 ^) ]$ D9 Q& c1 M, t3 N: p
3 |6 X# b. N  i5 B: V) c3 N  I

subplot(2,1,2);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

subplot(2,1,1);

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

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

xlabel('频率');

ylabel('幅度');

$ G9 v9 e* T2 w0 y* g
# @& _+ u7 @3 _3 }5 h% @' I

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;

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

1. q15_t testInput_radix2_q15_50hz[TEST_LENGTH_SAMPLES];
   / N: p) V0 R' B* l
2. q15_t testOutputQ15[TEST_LENGTH_SAMPLES/2];
   + g% s  d1 `7 o5 s: a. H+ _) k! b8 e
3. 
   1 a8 L6 C7 X. e& S
4. /*
   ; i+ |1 u: P8 y
5. *********************************************************************************************************
   & c5 Q) O1 O( A
6. *        函 数 名: arm_cfft_radix2_q15_app
   : f( o, e  H' @3 v+ c
7. *        功能说明: 调用函数arm_cfft_radix2_q15计算幅频。
   
8. *        形    参：无
   & r9 ~- s4 b; S9 }; _
9. *        返 回 值: 无
   
10. *********************************************************************************************************
    , @' g5 B8 a- x" C, p  M% Z( f+ M0 W5 @
11. */
    0 `( y2 V, G  m0 Q
12. void arm_cfft_radix2_q15_app(void)
    ; B" Y* {& u3 ~2 e2 G4 N# Y
13. {
    2 C( a. w2 Q9 w! d: u# C1 ]
14. uint16_t i,j;
    
15. arm_cfft_radix2_instance_q15 S;
    
16. fftSize = 1024;
    , s1 ], f% z5 {1 {9 Z# V. Z
17.     ifftFlag = 0;
    
18.     doBitReverse = 1;
    / m' q" S7 T5 Y7 A) N0 H! ?) R( L- w
19. /* 初始化结构S */
    
20. arm_cfft_radix2_init_q15(&S, fftSize, ifftFlag, doBitReverse);
    # U$ U: A4 o9 M' \' z8 ~) z
21. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    0 V# u8 ]7 X2 L2 ^
22. for(i=0; i<1024; i++)
    " [% y: L! v8 q4 }% ]) d, `
23. {
    ' ^+ u) l4 h# y% P/ @7 _% f
24. /* 虚部全部置0 */
    8 z  y! ~+ F: d' X" |$ c! D- L) |5 t
25. testInput_radix2_q15_50hz[i*2+1] = 0;
    
26. /* 51.2Hz正弦波，采样率1024Hz。
    
27.    arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
    
28.    32768 / 20 = 1638.4
    1 W3 a5 k- q7 d, t
29. */
    4 Y" m3 u# Y$ b! Y/ e
30. j = i % 20;
    
31.           testInput_radix2_q15_50hz[i*2] = arm_sin_q15(1638*j);
    7 W+ q, v0 j8 ?1 u5 R9 H
32.           printf("%d\r\n", testInput_radix2_q15_50hz[i*2]);
    0 q' f) ^' x. f+ t
33. }
    3 U( G8 e3 L. x1 y
34. /* 输出结果分割线 */
    
35. printf("**************************************************\r\n");
    
36. printf("**************************************************\r\n");
    
37. /* 计算CFFT */
    
38. arm_cfft_radix2_q15(&S, testInput_radix2_q15_50hz);
    # L  {* M- [7 B+ p4 H9 _3 N! F' g: c
39. /* 计算模值 */
    & N: G8 y8 |' [4 m+ C; H% E! y
40. arm_cmplx_mag_q15(testInput_radix2_q15_50hz, testOutputQ15, fftSize);
    8 X1 x" r) _; g7 R0 a
41. /* 串口打印求解的模值 */
    
42. for(i=0; i<1024; i++)
    
43. {
    $ D/ `2 F4 F% y7 z
44. printf("%d\r\n", testOutputQ15[i]);
    - F8 h0 u. e9 n( l! `! m
45. }
    ( e5 I5 J0 _" ?, i) I
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变换

6 f9 x6 x% i4 R7 Y4 {

subplot(2,1,2);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

subplot(2,1,1);

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

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

xlabel('频率');

ylabel('幅度');

  W) q2 z1 W1 \
$ ?' n! {) t3 U2 W: |% U, M/ \

Matlab运行结果如下：

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

' B: b. H/ a8 c# l9 H' @

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];
   
2. /*
   - Z1 K1 C8 e: Y2 d
3. *********************************************************************************************************
   
4. *        函 数 名: arm_cfft_radix4_q31_app
   
5. *        功能说明: 调用函数arm_cfft_radix4_q31_app计算幅频。
   2 Q7 @6 w7 L. L& a
6. *        形    参：无
   
7. *        返 回 值: 无
   
8. *********************************************************************************************************
   
9. */
   3 Q# w4 g% e# w+ _, p, _+ N& I
10. void arm_cfft_radix4_q31_app(void)
    7 G* m2 C7 L4 d* n$ E" O: }
11. {
    , J6 w: f3 c! _2 E1 E7 |
12. uint16_t i,j;
    
13. arm_cfft_radix4_instance_q31 S;
    
14. fftSize = 1024;
    
15.     ifftFlag = 0;
    
16.     doBitReverse = 1;
    % W+ B8 C( U7 D' i
17. /* 初始化结构S */
    8 [/ U5 e1 N6 v$ m/ u1 {# _) E
18. arm_cfft_radix4_init_q31(&S, fftSize, ifftFlag, doBitReverse);
    
19. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    
20. for(i=0; i<1024; i++)
    
21. {
    
22. testInput_radix4_q31_50hz[i*2+1] = 0;
    
23. /* 51.2Hz正弦波，采样率1024Hz。
    
24.    arm_sin_q31输入参数的范围0-2^31， 这里每20次为一个完整的正弦波，
    
25.    2^31 / 20 = 107374182.4
    
26. */
    
27. j = i % 20;
    
28.           testInput_radix4_q31_50hz[i*2] = arm_sin_q31(107374182*j);
    $ B1 i& h" p5 X7 m9 d
29.           printf("%d\r\n", testInput_radix4_q31_50hz[i*2]);
    
30. }
    ; U# m, [8 K# f1 U  `; ?* H3 p
31. /* 输出结果分割线 */
    
32. printf("**************************************************\r\n");
    
33. printf("**************************************************\r\n");
    
34. /* 计算CFFT */
    2 \* v, x) @: j$ I/ e; H2 d# O2 y( x5 R
35. arm_cfft_radix4_q31(&S, testInput_radix4_q31_50hz);
    
36. /* 计算模值 */
    
37. arm_cmplx_mag_q31(testInput_radix4_q31_50hz, testOutputQ31, fftSize);
    9 x; x9 ]  q. E
38. /* 串口打印求解的模值 */
    0 r* i6 H' J$ E" ]3 R/ f' R
39. for(i=0; i<1024; i++)
    ) M8 X9 \& M8 \* g
40. {
      N8 j3 |/ d" ~5 H' ~
41. printf("%d\r\n", testOutputQ31[i]);
    " {  u  J) z* e% t. P
42. }
    
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;           %真实的频率

, Y5 r2 C1 e9 A

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

' Z* I4 b/ x# b9 y1 d/ m4 a

subplot(2,1,2);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

- W/ \0 e+ c) w* _

subplot(2,1,1);

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

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

xlabel('频率');

ylabel('幅度');

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];
   
2. /*
   
3. *********************************************************************************************************
   / `7 y7 y! @- `" F; \/ R7 T& c
4. *        函 数 名: arm_cfft_radix4_q15_app
   
5. *        功能说明: 调用函数arm_cfft_radix4_q15计算幅频。
   # K: I9 |2 \/ i# v
6. *        形    参：无
   + v" i5 H$ B1 Z3 m! F
7. *        返 回 值: 无
   8 j, n+ h% Y! {$ m6 s
8. *********************************************************************************************************
   
9. */
   . J: k5 m3 _$ n
10. void arm_cfft_radix4_q15_app(void)
    1 k% P  R: ]8 G. ]
11. {
    
12. uint16_t i,j;
    $ z2 z6 m0 X6 u
13. arm_cfft_radix4_instance_q15 S;
    4 i6 J1 ]& Z6 U) G: y
14. fftSize = 1024;
    
15.     ifftFlag = 0;
    . B/ u7 e2 J8 J) j0 H% x1 o
16.     doBitReverse = 1;
    
17. /* 初始化结构S */
    ; g1 ?0 v8 g" u, E% U
18. arm_cfft_radix4_init_q15(&S, fftSize, ifftFlag, doBitReverse);
    
19. /* 按照实部，虚部，实部，虚部..... 的顺序存储数据 */
    
20. for(i=0; i<1024; i++)
    
21. {
    9 n- |; L- I1 ]
22. /* 虚部全部置0 */
    * g4 ~3 H2 F2 e% {
23. testInput_radix4_q15_50hz[i*2+1] = 0;
    
24. /* 51.2Hz正弦波，采样率1024Hz。
    6 j8 ]8 p5 d4 }& @: O, n. C8 q
25.    arm_sin_q15输入参数的范围[0, 32768)， 这里每20次为一个完整的正弦波，
    
26.    32768 / 20 = 1638.4
    
27. */
      u) c9 V. [$ C" \: f
28. j = i % 20;
    
29.           testInput_radix4_q15_50hz[i*2] = arm_sin_q15(1638*j);
    $ D! m+ U4 B1 M. i. q6 P7 k
30.           printf("%d\r\n", testInput_radix4_q15_50hz[i*2]);
    & d  K* o) _" k* [
31. }
    # }: a% F5 u( B" E# e' w
32. /* 输出结果分割线 */
    
33. printf("**************************************************\r\n");
    
34. printf("**************************************************\r\n");
    
35. /* 计算CFFT */
    4 [" f  r5 z! d' `6 ^" ]
36. arm_cfft_radix4_q15(&S, testInput_radix4_q15_50hz);
    " e3 ~2 u7 `% A9 z+ Z
37. /* 计算模值 */
    
38. arm_cmplx_mag_q15(testInput_radix4_q15_50hz, testOutputQ15, fftSize);
    # w+ ]! w2 c* F; f# c5 {
39. /* 串口打印求解的模值 */
    & X2 b; f' I1 N/ N2 B; S2 S4 R
40. for(i=0; i<1024; i++)
    
41. {
    
42. printf("%d\r\n", testOutputQ15[i]);
    ( _0 o3 q8 y2 J; O
43. }
    , s9 b% A  d5 [4 p% \, S
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;                %真实的频率

( G% P; ~* i5 f

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

0 ?5 Q+ d& V. M7 ?4 t

subplot(2,1,2);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

, V! {) e$ y6 \) L" l+ B1 @

subplot(2,1,1);

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

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

xlabel('频率');

ylabel('幅度');

, h$ @& s& O  G1 A8 ]

Matlab运行结果如下：

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

+ N9 y# D. }7 t# e

30.4 总结

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

! Q/ ^/ u. R; z0 _  W- t) I

标记下学习

mark一下！