特别说明：完整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 总结。

3 k0 q4 f* t4 z28.1 汇编FFT库说明
28.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库的移植就完成了。

! t+ g2 M% |. W- T6 ]6 o( y* g4 ?- p

) x& P5 n4 y( p; F
28.2 函数cr4_fft_1024_stm32的使用

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

1. 
   
2. uint32_t input[1024], output[1024], Mag[1024]; /* 输入，输出和幅值 */
   
3. /*
   
4. *********************************************************************************************************
   2 Z; L- V1 w% B' z0 |
5. *        函 数 名: PowerMag
   ' s4 ?. q2 O: b& Q! Q# z  _  K
6. *        功能说明: 求模值
   0 k1 A9 ~/ p( y  M3 K3 t/ Q
7. *        形    参：_usFFTPoints  FFT点数
   
8. *        返 回 值: 无
   
9. *********************************************************************************************************
   
10. */
    
11. void PowerMag(uint16_t _usFFTPoints)
    7 p4 M( l/ _# W  ^. [3 d
12. {
    - s9 @7 x2 |, t" g1 d2 K! @% [
13. int16_t lX,lY;
    
14. uint16_t i;
    
15. float32_t mag;
    
16. 
    
17. /* 计算幅值 */
    
18. for (i=0; i < _usFFTPoints; i++)
    
19. {
    6 D" i- V, Z! X' m" @, B
20.           lX= (output[i]<<16)>>16;          /* 实部*/
    
21. lY= (output[i]>> 16);             /* 虚部 */   
    
22. mag = __sqrtf(lX*lX+ lY*lY);      /* 求模 */
    
23. Mag[i]= mag*2;                    /* 求模后乘以2才是实际模值，直流分量不需要乘2 */
    
24. }
    
25.      
    
26. /* 由于上面多乘了2，所以这里直流分量要除以2 */
    
27. Mag[0] = Mag[0]>>1;
    3 m* U4 F' Y! Q- I: v
28. }
    7 g3 S' t6 {8 R( T$ o
29. 
    ; i* D5 J7 L9 B6 p& _4 S
30. /*
    8 x: ]& {) c7 {0 O$ r* m
31. *********************************************************************************************************
    ; e( U. t2 D" Z! l9 ^6 @
32. *        函 数 名: DSP_FFT1024
    
33. *        功能说明: 1024点FFT实现
    ) ~% G: B4 t" z& s
34. *        形    参：无
    
35. *        返 回 值: 无
    , Z2 J) y* j9 \9 f# q  v
36. *********************************************************************************************************
    6 U/ F2 ], n* C4 Z+ I5 U( t- X
37. */
    " u1 K  C) \: }8 B$ K8 ]
38. void DSP_FFT1024(void)
    8 }( K# w3 s6 `! N2 A
39. {
    ) Y( p* E8 t: P4 N) o
40. uint16_t i;
    
41. /* 获得1024个采样点 */
    8 n) T' Y. V% [' p0 d
42. for (i = 0; i < 1024; i++)
    4 u- H6 v! |$ L' L+ Z
43. {
    0 e5 H- P  L9 G( l
44. input[i] = 0;
    , D$ W2 n3 l! W5 z; m! c! @% f
45. /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率1KHz */
    " S- a' B7 ^# v
46. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000) + 512*sin(2*3.1415926f*20*i/1000) ;
    
47. }
    
48. /* 计算1024点FFT
    4 z- c' |0 H! \  ?1 P0 V+ \# v
49.    output：输出结果，高16位是虚部，低16位是实部。
    . Y/ X* P9 q4 M" p
50.    input ：输入数据，高16位是虚部，低16位是实部。
    0 x6 [& E+ ?. F5 o/ i! Y# Y
51.    第三个参数必须是1024。
    
52. */
    
53. cr4_fft_1024_stm32(output, input, 1024);
    
54. /* 求幅值 */
    
55. PowerMag(1024);
    - r/ z5 B6 J3 k/ o" }1 w' O1 P8 w( I- h
56. /* 打印输出结果 */
    
57. for (i = 0; i < 1024; i++)
    
58. {
    " V, _( k& }: B; X" X
59. printf("%d\r\n", Mag[i]);
    
60. }
    
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;          %真实的频率

& x$ `- g$ ]( I* p+ S
. r9 X9 ^/ o8 G0 Y5 @& I

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

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

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

" `* e7 y) \: P1 K9 }4 j# ~

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计算的结果频谱泄露稍严重些（具体情况要具体分析）。

6 {. `) \, V# i7 [; D! V1 F28.3 函数cr4_fft_256_stm32的使用

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

1. /*
   
2. *********************************************************************************************************
   7 c1 |1 P3 W+ @  @. p8 U
3. *        函 数 名: DSP_FFT256
   
4. *        功能说明: 256点FFT实现
   7 K/ Y* b, ~' Z7 X* m" o4 h
5. *        形    参：无
   
6. *        返 回 值: 无
   ; n( B4 |8 @  j" Q  O' L
7. *********************************************************************************************************
   
8. */
   6 K  I0 d5 u. c- b1 @3 u
9. void DSP_FFT256(void)
   ' x' Y% ], r8 u. \. ~
10. {
    
11. uint16_t i;
    
12. /* 获得256个采样点 */
    ' [0 X: G9 t9 v1 T
13. for (i = 0; i < 256; i++)
    
14. {
    
15. input[i] = 0;
    
16. /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率200Hz */
    4 {( D# m% L6 [6 }  r$ u
17. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/200) + 512*sin(2*3.1415926f*20*i/200) ;
    
18. }
    
19. /* 计算256点FFT
    
20.    output：输出结果，高16位是虚部，低16位是实部。
    
21.    input ：输入数据，高16位是虚部，低16位是实部。
    
22.    第三个参数必须是1024。
    
23. */
    
24. cr4_fft_256_stm32(output, input, 256);
    
25. /* 求幅值 */
    
26. PowerMag(256);
    
27. /* 打印输出结果 */
    9 y( C) d- x0 d) g3 h8 U; ^# D: G
28. for (i = 0; i < 256; i++)
    " j. d5 R5 w8 c' \% K6 e% P& X$ R
29. {
    & w- z3 _' ]! v  c! O2 @
30. printf("%d\r\n", Mag[i]);
    . D* p0 N$ x( Z6 e1 [
31. }
    % a" o3 ~5 V% N& T7 I3 r7 ^1 {
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;          %真实的频率

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

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

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

+ i/ }, T/ q1 [
4 C% Q4 m+ E! d: N/ L2 E

subplot(2,1,1);

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

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

  P4 r  S: q7 e/ L
9 b* T* f; P1 l7 W& f

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的效果都比较一般。

4 C3 E2 B: r% J% S/ u

28.4 函数cr4_fft_64_stm32的使用

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

1. /*
   
2. *********************************************************************************************************
   / s, h, o6 h  m" F/ j$ i
3. *        函 数 名: DSP_FFT64
   ) d9 f' i' H9 {% G/ B- P3 d% f
4. *        功能说明: 64点FFT实现
   6 Q; J; _$ E# l- H. Z9 u
5. *        形    参：无
   % V3 T7 L% V2 ?+ U$ i, e
6. *        返 回 值: 无
   
7. *********************************************************************************************************
   3 Z# [5 x6 f* J4 \# O5 X8 O& g2 e
8. */
   ; Y( U  ^+ M# X  P8 h
9. void DSP_FFT64(void)
   2 U7 r! [! o( K  ?, R; f
10. {
    / V2 K2 t& V, }3 X
11. uint16_t i;
    
12. /* 获得64个采样点 */
    
13. for (i = 0; i < 64; i++)
    # D# L! O. ~8 w( D7 H7 w
14. {
    1 R6 Y' @1 R' ?2 t5 b1 V. m( b
15. input[i] = 0;
    
16. /* 波形是由直流分量，5Hz正弦波和10Hz正弦波组成，波形采样率60Hz */
    
17. input[i] = 1024 + 1024*sin(2*3.1415926f*5*i/60) + 512*sin(2*3.1415926f*10*i/60) ;
    
18. }
    , j: A0 W; ?0 f7 Z* j- f" J
19. /* 计算64点FFT
    # Q* B" T4 y- y4 y
20.    output：输出结果，高16位是虚部，低16位是实部。
    ( j% N- t. l' i/ I
21.    input ：输入数据，高16位是虚部，低16位是实部。
    
22.    第三个参数必须是1024。
    2 w) m& q9 q2 @* V
23. */
    & \; n" V! N! Q% H( Q
24. cr4_fft_64_stm32(output, input, 64);
    
25. /* 求幅值 */
    % q4 t3 L% {( N8 k) _7 j) Y& T
26. PowerMag(64);
    : A; }5 _2 O7 |$ e+ {4 K3 v
27. /* 打印输出结果 */
    
28. for (i = 0; i < 64; i++)
    # t. F1 ~' c# s" G% ?( c
29. {
    ) c% [1 ]4 b  x8 ~3 {  t
30. printf("%d\r\n", Mag[i]);
    6 k$ O3 i+ F( m4 Y7 @& Q2 i3 E" h
31. }
    
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;          % 真实的频率

% Z0 s* t- [5 N" X+ o

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

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

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

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. /*
   7 h' f# k, F. {0 l, Q, x
2. *********************************************************************************************************
   ! U& P7 o4 J+ }' y$ S
3. *        函 数 名: Power_Phase_Radians
   
4. *        功能说明: 求相位
   
5. *        形    参：_usFFTPoints  FFT点数
   
6. *        返 回 值: 无
   
7. *********************************************************************************************************
   6 u% Q& }. Q9 `; B) v: O
8. */
   6 f, t. V; d0 C! s
9. void Power_Phase_Radians(uint16_t _usFFTPoints)
   4 J6 n* Q/ a# T
10. {
    
11. int16_t lX, lY;
    
12. uint16_t i;
    1 `6 V. b: v: \! S! i( h0 M$ {8 ?
13. float32_t phase;
    # ]! }7 H; O3 z' f5 U
14. for (i=0; i <_usFFTPoints; i++)
    
15. {
    
16. lX= (output[i]<<16)>>16;  /* 实部 */
    + k$ A3 r) V; H; Z
17. lY= (output[i] >> 16);    /* 虚部 */
    
18.           phase = atan2(lY, lX);    /* atan2求解的结果范围是(-pi, pi], 弧度制 */
    
19.    
    
20. Phase[i] = phase* 180.0f/3.1415926f;  /* 将求解的结果由弧度转换为角度 */
    & x# A1 h; T: t+ ^5 D4 L1 Q1 t( s- j
21. }
    , O  U1 d3 O4 K  l4 r
22. }
    
23. 
    
24. /*
    ) E6 e: s( a# N. l
25. *********************************************************************************************************
    ! e8 A% h" v/ W
26. *        函 数 名: DSP_FFTPhase
    
27. *        功能说明: 1024点FFT的相位求解
    / _4 s! I1 [; }) m: V
28. *        形    参：无
    ) |5 u% O9 N, G0 V4 y4 D# N# e' A+ p
29. *        返 回 值: 无
    
30. *********************************************************************************************************
    
31. */
    
32. void DSP_FFTPhase(void)
    
33. {
    
34. uint16_t i;
    
35. /* 获得1024个采样点 */
    ! ~8 T1 V& A/ X! h: h/ n/ a
36. for (i = 0; i < 1024; i++)
    5 w7 S3 ?2 f3 D6 ?7 ]$ D7 v
37. {
    
38. input[i] = 0;
    
39. /* 波形是由直流分量，50Hz正弦波组成，波形采样率1KHz */
    " Z; K, n" R# k1 H: y7 z& ]
40. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000 + 3.1415926f/3);
    & i. O0 E0 G; t9 [4 X# _. [
41. }
    5 }9 ]7 @0 q3 k- Z' O( [
42. /* 计算1024点FFT
    + E% n. r% \* G! o! j) l7 I
43.    output：输出结果，高16位是虚部，低16位是实部。
    
44.    input ：输入数据，高16位是虚部，低16位是实部。
    
45.    第三个参数必须是1024。
    
46. */
    0 M* i0 \# ^3 r" i8 d
47. cr4_fft_1024_stm32(output, input, 1024);
    # C( |( @& X6 e
48. /* 求相频 */
    
49. Power_Phase_Radians(1024);
    
50. /* 打印输出结果 */
    
51. for (i = 0; i < 1024; i++)
    
52. {
    0 u7 l/ c6 z/ U8 V/ a
53. printf("%f\r\n", Phase[i]);
    ! H; _# e; x' \/ R; b# D. Y
54. }
      e' v$ u$ ~: h
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;          %真实的频率

' V, D9 V. M7 k" Q8 q  O' a* v
6 y: z' P* y7 `5 _0 [

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

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

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

6 e8 [4 y4 b: O' |
  P8 @. ^9 F- d5 B. e! \  S

subplot(2,1,1);

realvalue = real(y);

imagvalue = imag(y);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

5 @! \' i2 h3 e1 L3 U# d# `

subplot(2,1,2);

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

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

* z3 {- |5 ]+ L2 B

Matlab运行结果如下：

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

9 O3 N% K! E3 e6 q; w2 A

28.6总结
本章节主要讲解了汇编FFT的1024点，256点和64点使用方法，有兴趣的可以深入了解汇编代码的实现。

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

标记学习

学习一下。。。

学习一下！！

) ~1 F1 j' p2 q8 F5 {学习一下。。。收藏了

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