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

第28章 ST官方汇编FFT库应用

    本章主要讲解ST官方汇编FFT库的应用，包括1024点，256点和64点FFT的实现。

    28.1 汇编FFT库说明

    28.2 函数cr4_fft_1024_stm32的使用

    28.3 函数cr4_fft_256_stm32的使用

    28.4 函数cr4_fft_64_stm32的使用

    28.5 汇编FFT的相频响应求解

    28.6 总结。

28.1 汇编FFT库说明
7 }$ Q, N* f4 M+ u9 G2 N28.1.1 描述

    这个汇编的FFT库是来自STM32F10x DSP library，由于是汇编实现的，而且是基4算法，所以实现FFT在速度上比较快。

    如果x[N]是采样信号的话，使用FFT时必须满足如下两条：

l N得满足 （n =1，2, 3…..），也就是以4为基数。

l 采样信号必须是32位数据，高16位存实部，低16位存虚部（这个是针对大端模式），小端模式是高位存虚部，低位存虚部。一般常用的是小端模式。

汇编FFT的实现主要包括以下三个函数：

    1. cr4_fft_64_stm32 ：实现64点FFT。

    2. cr4_fft_256_stm32 ：实现256点FFT。

    3. cr4_fft_1024_stm32 ： 实现1024点FFT。

28.1.2 汇编库的移植

    这个汇编库的移植比较简单，用户需要从网上搜索STM32F10x_DSP_Lib_V2.0.0（官网没有找到这个软件包，所以需要用户在百度或者谷歌上搜索下）。下载后解压，在路径：STMicroelectronics\STM32F10x_DSP_Lib_V2.0.0\Libraries\STM32F10x_DSP_Lib\src\asm\arm下会看到如下文件：

    上面的是源文件，使用源文件还需要添加相应的头文件，头文件在路径：STMicroelectronics\STM32F10x_DSP_Lib_V2.0.0\Libraries\STM32F10x_DSP_Lib\inc下，文件如下所示：

    具备这几个文件就可以移植使用了，移植非常简单，把源文件的三个FFT库和两个头文件添加上即可，添加后效果如下（记得添加头文件的路径）：

    相应文件添加后还有最重要一条，要把stm32_dsp.h文件中的STM32F1头文件换成STM32F4的头文件：

经过上面的操作，汇编FFT库的移植就完成了。

% M6 w; O, r: z1 |
; z6 l4 d) U+ o" F+ B28.2 函数cr4_fft_1024_stm32的使用

    cr4_fft_1024_stm32用于实现1024点数据的FFT计算。下面通过在开发板上运行这个函数并计算幅频相应，然后再与Matlab计算的结果做对比。

1. 
   " t% @/ k  q8 v. K* I1 `/ k
2. uint32_t input[1024], output[1024], Mag[1024]; /* 输入，输出和幅值 */
   / p+ B4 K6 c: _& h
3. /*
   & _' T1 M" K, ?: C0 s  F
4. *********************************************************************************************************
   
5. *        函 数 名: PowerMag
   3 `# f2 ~$ ~1 x) o5 D  t
6. *        功能说明: 求模值
   " l$ ]6 Y0 ~* G. X6 Q
7. *        形    参：_usFFTPoints  FFT点数
     D: Q& A# r2 o- H
8. *        返 回 值: 无
   
9. *********************************************************************************************************
   
10. */
    
11. void PowerMag(uint16_t _usFFTPoints)
    
12. {
    
13. int16_t lX,lY;
    
14. uint16_t i;
    
15. float32_t mag;
    
16. 
    
17. /* 计算幅值 */
    1 g5 C2 S7 R6 Y8 {3 t# |
18. for (i=0; i < _usFFTPoints; i++)
    & j% Z3 i2 A5 o. \$ L% D8 E
19. {
    / h$ }9 M, @/ c* F& _
20.           lX= (output[i]<<16)>>16;          /* 实部*/
    
21. lY= (output[i]>> 16);             /* 虚部 */   
    
22. mag = __sqrtf(lX*lX+ lY*lY);      /* 求模 */
    ; [5 E4 C# R$ N" r# Y4 q3 Z# \/ ~
23. Mag[i]= mag*2;                    /* 求模后乘以2才是实际模值，直流分量不需要乘2 */
    
24. }
    
25.      
    ; y. o1 h0 q+ q/ A& O: ~
26. /* 由于上面多乘了2，所以这里直流分量要除以2 */
    3 q4 _8 B8 ?" u6 I3 O8 h3 O/ d
27. Mag[0] = Mag[0]>>1;
    - f( w. m) u6 h( M+ w. B2 p
28. }
    ; I& `" P. p# ?' I
29. 
    ( a" O; q$ s- v
30. /*
    0 w/ N& f. u: M6 g) c- d  j
31. *********************************************************************************************************
    
32. *        函 数 名: DSP_FFT1024
    
33. *        功能说明: 1024点FFT实现
    : g9 d: ~5 e( A1 o$ _
34. *        形    参：无
    
35. *        返 回 值: 无
    
36. *********************************************************************************************************
    ) B  F. @+ M+ Y+ p7 O+ m0 i$ \; Z- c
37. */
    , ~4 g" {9 E) l1 r1 E
38. void DSP_FFT1024(void)
    0 V( T# J+ J+ N. a
39. {
    ( ?/ p0 h) N8 n* Z0 q# A
40. uint16_t i;
    
41. /* 获得1024个采样点 */
    
42. for (i = 0; i < 1024; i++)
    3 p6 Q- O: Q% ], K8 g
43. {
    
44. input[i] = 0;
    * q! T4 e5 ^9 P, V
45. /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率1KHz */
    
46. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000) + 512*sin(2*3.1415926f*20*i/1000) ;
    
47. }
    
48. /* 计算1024点FFT
    - g: R6 `1 i' R9 ^* E. s& e- E9 V# b
49.    output：输出结果，高16位是虚部，低16位是实部。
    + W8 B5 h6 A5 Q& C# Q
50.    input ：输入数据，高16位是虚部，低16位是实部。
    
51.    第三个参数必须是1024。
      [( r  F% a# r. C2 d
52. */
    
53. cr4_fft_1024_stm32(output, input, 1024);
    
54. /* 求幅值 */
    
55. PowerMag(1024);
    
56. /* 打印输出结果 */
    , M0 g* Y! u! h* {% d
57. for (i = 0; i < 1024; i++)
    
58. {
    
59. printf("%d\r\n", Mag[i]);
    7 Q# z9 h% D* O0 X5 ?' ~
60. }
    * q( J/ a. ^8 }* x
61. }

复制代码

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

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

Fs = 1000;                  % 采样率

N  = 1024;           % 采样点数

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

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

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

/ D2 @- [  [1 S7 V5 ?  D

%波形是由直流分量，50Hz正弦波和20Hz正弦波组成

x = 1024 + 1024*sin(2*pi*50*t) + 512*sin(2*pi*20*t)  ;  

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

7 t7 ?3 F' }1 B! A1 v  q

subplot(2,1,1);

Mag = abs(y)*2/N;         %求FFT转换结果的模值

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

subplot(2,1,2);

plot(f, sampledata);   %绘制STM32计算的幅频相应

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

    运行Matlab后的输出结果如下：

    从上面的对比结果中可以看出，Matlb和函数cr4_fft_1024_stm32计算的结果基本是一直的，但是函数cr4_fft_1024_stm32计算的结果频谱泄露稍严重些（具体情况要具体分析）。

+ b$ [7 f( e- J3 E& \, F
8 C5 X, G3 P% M3 h28.3 函数cr4_fft_256_stm32的使用

    cr4_fft_256_stm32和cr4_fft_1024_stm32的用法是一样的，下面通过一个实例进行说明：

1. /*
   8 ~. u6 L1 ^4 y9 z+ L3 F3 i( ]
2. *********************************************************************************************************
   
3. *        函 数 名: DSP_FFT256
   ( r9 `2 B# f/ N0 d; R0 n% C  c
4. *        功能说明: 256点FFT实现
   
5. *        形    参：无
   
6. *        返 回 值: 无
   
7. *********************************************************************************************************
   , |! _/ G, p7 G# |8 U  x6 y
8. */
   
9. void DSP_FFT256(void)
   0 S2 a5 g' H7 l2 }6 ~' @$ }
10. {
    
11. uint16_t i;
    
12. /* 获得256个采样点 */
    
13. for (i = 0; i < 256; i++)
    % \! R! y; t0 v% Q0 F
14. {
    
15. input[i] = 0;
    ; ?5 B# i1 `; W1 G) y
16. /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率200Hz */
    
17. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/200) + 512*sin(2*3.1415926f*20*i/200) ;
    
18. }
    
19. /* 计算256点FFT
    ' r5 o2 t$ _, d5 ^
20.    output：输出结果，高16位是虚部，低16位是实部。
    2 h1 P& q: i. ?1 Y) O- {+ l
21.    input ：输入数据，高16位是虚部，低16位是实部。
    0 ]  B( C8 u( Z$ p( `
22.    第三个参数必须是1024。
    
23. */
    ) r. _& j3 A' i& ^: x8 [
24. cr4_fft_256_stm32(output, input, 256);
      d6 F: ]7 }; ~( F7 t" h) c
25. /* 求幅值 */
    
26. PowerMag(256);
    
27. /* 打印输出结果 */
    , `" G7 w: X" W
28. for (i = 0; i < 256; i++)
    6 m+ x2 f  \# _2 B- z. v" ^; x
29. {
    ) _' `0 n/ E8 `
30. printf("%d\r\n", Mag[i]);
    $ B4 M7 p- y6 V3 N4 w: Y4 ]( e
31. }
    ! q+ N6 S0 ?5 f( d
32. }

复制代码

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

    对比前需要先将串口打印出的数据加载到Matlab中，并给这个数组起名sampledata， Matlab中运行的代码如下：

Fs = 200;              % 采样率

N  = 256;             % 采样点数

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

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

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

. j4 a4 y' q% _4 l7 O( t6 e8 r

%波形是由直流分量，50Hz正弦波和20Hz正弦波组成

x = 1024 + 1024*sin(2*pi*50*t) + 512*sin(2*pi*20*t)  ;  

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

5 r" J$ f4 U6 r& p& {

subplot(2,1,1);

Mag = abs(y)*2/N;    %求FFT转换结果的模值

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

4 g* P7 \5 t2 z

subplot(2,1,2);

plot(f, sampledata);   %绘制STM32计算的幅频相应

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

    运行Matlab后的输出结果如下：

    从上面的对比结果中可以看出，Matlb和函数cr4_fft_1024_stm32计算的结果基本是一直的，但是函数cr4_fft_1024_stm32计算正弦波中心频率更好些，左右两边的频谱泄露较少，但由于计算的点数较少，Matlb和函数cr4_fft_1024_stm32的效果都比较一般。

28.4 函数cr4_fft_64_stm32的使用

    cr4_fft_256_stm32和cr4_fft_1024_stm32的用法是一样的，下面通过一个实例进行说明：

1. /*
   * B$ R0 @# S1 u. [/ l1 }0 x
2. *********************************************************************************************************
   
3. *        函 数 名: DSP_FFT64
   
4. *        功能说明: 64点FFT实现
   
5. *        形    参：无
   
6. *        返 回 值: 无
   1 g5 M1 O& R1 {- q- b
7. *********************************************************************************************************
   
8. */
   
9. void DSP_FFT64(void)
   / N0 d3 b) F+ a4 s
10. {
    
11. uint16_t i;
    
12. /* 获得64个采样点 */
    
13. for (i = 0; i < 64; i++)
    
14. {
    
15. input[i] = 0;
    
16. /* 波形是由直流分量，5Hz正弦波和10Hz正弦波组成，波形采样率60Hz */
    0 J& E& K; O6 ?9 A: k
17. input[i] = 1024 + 1024*sin(2*3.1415926f*5*i/60) + 512*sin(2*3.1415926f*10*i/60) ;
    , c) s, Z! a- s$ V% p8 V
18. }
    ) a# _7 Q) w  C( k$ e
19. /* 计算64点FFT
    / D' Z( Z% u" ?% o
20.    output：输出结果，高16位是虚部，低16位是实部。
    
21.    input ：输入数据，高16位是虚部，低16位是实部。
    
22.    第三个参数必须是1024。
    , ?0 L5 A* j. y  x1 B( l
23. */
    
24. cr4_fft_64_stm32(output, input, 64);
    + R1 j+ t9 v, f3 H9 A5 ]) c
25. /* 求幅值 */
    
26. PowerMag(64);
    9 W) T3 c6 d( [
27. /* 打印输出结果 */
    3 t- E& x' y7 q
28. for (i = 0; i < 64; i++)
    ! g6 B0 m6 Z4 h7 R
29. {
    + I; u, Y1 Z) s7 {' O* m6 p
30. printf("%d\r\n", Mag[i]);
    ( y1 K9 o1 a8 T# r. D: b
31. }
    , {0 R6 f; g) J
32. }

复制代码

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

    对比前需要先将串口打印出的数据加载到Matlab中，并给这个数组起名sampledata，Matlab中运行的代码如下：

Fs = 60;               % 采样率

N  = 64;              % 采样点数

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

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

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

%波形是由直流分量，5Hz正弦波和10Hz正弦波组成

x = 1024 + 1024*sin(2*pi*5*t) + 512*sin(2*pi*10*t)  ;  

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

  J+ v# }4 |# E1 L: A. ?. f" w

subplot(2,1,1);

Mag = abs(y)*2/N;    %求FFT转换结果的模值

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

subplot(2,1,2);

plot(f, sampledata);   %绘制STM32计算的幅频相应

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

    运行Matlab后的输出结果如下：

    从上面的对比结果中可以看出，Matlb和函数cr4_fft_64_stm32计算的结果基本是一直的，但是计算的效果都比较差，主要是因为采样点数太少。

28.5 汇编FFT的相频响应求解

    特别说明：这里相频响应求解暂时有问题（与Matlab求解的结果不一致），但是也把这个求解过程

列出来作为大家的参考。

    相频响应代码求解如下：

1. /*
   " {# s& i- n9 z* q7 ]( A) M& R
2. *********************************************************************************************************
     u# [' l" I1 m
3. *        函 数 名: Power_Phase_Radians
   5 y7 s- Y' i3 b: A8 z- K) J" n- u
4. *        功能说明: 求相位
   
5. *        形    参：_usFFTPoints  FFT点数
   
6. *        返 回 值: 无
   ( J( U( q1 X* T3 Q: j6 v
7. *********************************************************************************************************
   
8. */
     B$ _, Z- @7 v* H
9. void Power_Phase_Radians(uint16_t _usFFTPoints)
   . F8 X1 R2 \$ `: i
10. {
    
11. int16_t lX, lY;
    
12. uint16_t i;
    ) |1 s# w* X0 D8 E
13. float32_t phase;
    
14. for (i=0; i <_usFFTPoints; i++)
    ; S6 p; g! A* p8 u2 s: Q
15. {
    " |7 `8 H# N/ _% w$ _) A) f
16. lX= (output[i]<<16)>>16;  /* 实部 */
    9 \$ s4 d, i, S6 h* `2 _9 G
17. lY= (output[i] >> 16);    /* 虚部 */
    
18.           phase = atan2(lY, lX);    /* atan2求解的结果范围是(-pi, pi], 弧度制 */
    4 c: m. `  b( X0 S  ]% c% _
19.    
    
20. Phase[i] = phase* 180.0f/3.1415926f;  /* 将求解的结果由弧度转换为角度 */
    3 |* a* z- d& H, J0 u, G. e+ p) z
21. }
    4 ~$ `' v# E, {, c* K( E
22. }
    6 Q3 Q' `4 j; {. c. o8 I
23. 
    * c( w& ?# u# T
24. /*
    
25. *********************************************************************************************************
    
26. *        函 数 名: DSP_FFTPhase
    
27. *        功能说明: 1024点FFT的相位求解
    
28. *        形    参：无
    
29. *        返 回 值: 无
      ]! x: b1 W6 Y* P
30. *********************************************************************************************************
    3 i" g5 U7 ^" R' r: V1 h
31. */
    , s+ ?! _7 N* b7 {# Y
32. void DSP_FFTPhase(void)
    
33. {
    
34. uint16_t i;
    ; L7 y; J) E+ q4 ]' [
35. /* 获得1024个采样点 */
    
36. for (i = 0; i < 1024; i++)
    
37. {
    
38. input[i] = 0;
    
39. /* 波形是由直流分量，50Hz正弦波组成，波形采样率1KHz */
    ) R5 b* a! m& x
40. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000 + 3.1415926f/3);
    # ?& Z. V$ `4 e5 Z- w- N* G$ ]
41. }
    
42. /* 计算1024点FFT
    
43.    output：输出结果，高16位是虚部，低16位是实部。
    5 A9 K: d% e* W, V
44.    input ：输入数据，高16位是虚部，低16位是实部。
    7 ?8 A4 V, G8 C6 {0 w" l
45.    第三个参数必须是1024。
    2 Z! c* c. B' l6 E& K* I4 m4 i9 b6 K
46. */
    
47. cr4_fft_1024_stm32(output, input, 1024);
    . c, T" p1 k. |6 _' U2 Q
48. /* 求相频 */
    
49. Power_Phase_Radians(1024);
    5 B2 W# i4 f/ P8 F5 @: g
50. /* 打印输出结果 */
    ( q, k5 K3 u$ Z5 y7 C4 E: S6 _
51. for (i = 0; i < 1024; i++)
    ) V% m- `$ b! K5 B% z& W5 n8 H
52. {
    0 q+ F5 S0 u5 C1 _) V; T  }5 C
53. printf("%f\r\n", Phase[i]);
    % A- j4 w! Q( |! F7 L- y$ n
54. }
    
55. }

复制代码

运行函数DSP_FFTPhase可以通过串口打印出计算的相位，下面我们就通过Matlab计算的相位跟cr4_fft_1024_stm32计算的相位做对比。

    对比前需要先将串口打印出的数据加载到Matlab中，并给这个数组起名sampledata，Matlab中运行的代码如下：

Fs = 1000;             % 采样率

N  = 1024;           % 采样点数

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

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

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

%波形是由直流分量，50Hz正弦波正弦波组成

x = 1024 + 1024*sin(2*pi*50*t + pi/3)   ;  

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

# N) L  \! d$ I' J/ r' L7 t1 r: h
2 ~- K  H* z$ z

subplot(2,1,1);

realvalue = real(y);

imagvalue = imag(y);

plot(f, atan2(imagvalue, realvalue)*180/pi);

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

% \: I0 C7 n* d

subplot(2,1,2);

plot(f, sampledata);   %绘制STM32计算的幅频相应

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

Matlab运行结果如下：

从上面的运行结果来了，两者所得结果差距较大。有待后面进一步研究。

" M2 P% Q* s/ o+ i( K- w

28.6总结
7 G* p: u8 Z# e6 M7 @本章节主要讲解了汇编FFT的1024点，256点和64点使用方法，有兴趣的可以深入了解汇编代码的实现。

帮顶啊！！！！！！！！！！！！！！！！！！！！！！！！！！！

标记学习

学习一下。。。

学习一下！！

% m0 i- e- r0 |% ]) d学习一下。。。收藏了

这个是移植到M4平台的么,M7平台是否也可以这么用?