特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：链接
0 l( t8 F% e! R4 ^/ I$ R! N

第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库说明
3 }; q$ H+ r) J9 b: S8 q28.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库的移植就完成了。

! _9 x, @) \0 q6 k0 m5 ~

28.2 函数cr4_fft_1024_stm32的使用

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

1. 
   
2. uint32_t input[1024], output[1024], Mag[1024]; /* 输入，输出和幅值 */
   ) R' N& D$ ^% @- [  K6 ^
3. /*
   9 L1 n+ L1 [* }& l2 V  _
4. *********************************************************************************************************
   8 _+ D& P- p% |
5. *        函 数 名: PowerMag
   . k. s. {/ K- U; D
6. *        功能说明: 求模值
   
7. *        形    参：_usFFTPoints  FFT点数
   
8. *        返 回 值: 无
   
9. *********************************************************************************************************
     ^- f/ Q$ J# y1 z- m: _9 t
10. */
    
11. void PowerMag(uint16_t _usFFTPoints)
    , R. s6 I, c9 p3 `$ N
12. {
    $ W* ]& f) i% R' |- A
13. int16_t lX,lY;
    
14. uint16_t i;
    $ \2 w; r9 F$ M  m8 u0 ^& Y! E
15. float32_t mag;
    2 n# a: V6 x% i3 H- Y- T
16. 
    ; a/ O3 `1 q/ r7 x- |
17. /* 计算幅值 */
    
18. for (i=0; i < _usFFTPoints; i++)
    ' O) k5 ?2 j" M! l: i
19. {
    
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 */
    * T' M: D% S5 h! d
27. Mag[0] = Mag[0]>>1;
    6 {% e  R+ K) }2 t) ^7 N+ r
28. }
    
29. 
    
30. /*
    
31. *********************************************************************************************************
    
32. *        函 数 名: DSP_FFT1024
    ( r( R4 u3 V% i6 {
33. *        功能说明: 1024点FFT实现
    
34. *        形    参：无
    7 q1 U/ U- M7 O" I1 H6 M
35. *        返 回 值: 无
    . V5 g6 c& x* ^4 D
36. *********************************************************************************************************
    / j, G# w$ A; K  _" `5 w- f
37. */
    
38. void DSP_FFT1024(void)
    
39. {
    & a1 t9 O" k" e- V
40. uint16_t i;
    
41. /* 获得1024个采样点 */
      m0 Q, t7 z3 i) s6 f+ \- V
42. for (i = 0; i < 1024; i++)
    , O; X1 C1 r8 Y. K/ b; f0 M) z& f  n0 o
43. {
    1 x9 ?& i4 ^# s& I5 A4 o: X8 i3 u1 K
44. input[i] = 0;
    % `) i& w. J4 f( E: t- c* q
45. /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率1KHz */
    
46. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000) + 512*sin(2*3.1415926f*20*i/1000) ;
    ; c+ ?$ p- z9 o1 m( n
47. }
    
48. /* 计算1024点FFT
    0 S, z/ N% d# q
49.    output：输出结果，高16位是虚部，低16位是实部。
    7 B) F$ h5 K; M+ R; U( w) {' C
50.    input ：输入数据，高16位是虚部，低16位是实部。
    7 c$ h' f- U) m9 x. t
51.    第三个参数必须是1024。
    $ d1 M9 b1 U; i" Z2 _$ _& Z& h3 W$ ?7 e
52. */
    7 t" o/ b- j+ k; l
53. cr4_fft_1024_stm32(output, input, 1024);
    
54. /* 求幅值 */
    ! P; T% h; w8 u( i  Q' D
55. PowerMag(1024);
    & K. s* ?8 l) b4 \( ~8 t
56. /* 打印输出结果 */
    - T! `, a# {! _3 X* ^. E% A4 v
57. for (i = 0; i < 1024; i++)
    
58. {
    
59. printf("%d\r\n", Mag[i]);
    
60. }
    - o- B, c* G  M; W$ \# o0 Z( M
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;          %真实的频率

! F. Y8 x0 X- k  U7 q3 E* D3 f

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

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

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

* U, ]7 ^0 z" p3 R4 ^+ ~

subplot(2,1,1);

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

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

) o5 b0 q1 a  }  I5 t, a
; D. g9 P4 ~! ]: s8 s# D1 s

subplot(2,1,2);

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

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

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

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

7 A2 i7 s! ^  r9 t, n5 }& a. h

! Y8 `" U# w  T$ o7 V7 x28.3 函数cr4_fft_256_stm32的使用

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

1. /*
   2 [' X$ ]% ]9 G- s5 V1 `
2. *********************************************************************************************************
   4 w' i) l2 d0 R& E
3. *        函 数 名: DSP_FFT256
   " B4 @2 N& q, S2 k
4. *        功能说明: 256点FFT实现
   6 Q$ W4 ?  `# ?) J+ o# u6 x
5. *        形    参：无
   
6. *        返 回 值: 无
   + g; M: P3 X* b+ Y  u0 y
7. *********************************************************************************************************
   0 H  P6 }# u5 Z7 Q5 w
8. */
   
9. void DSP_FFT256(void)
   
10. {
    
11. uint16_t i;
    
12. /* 获得256个采样点 */
    + j$ x5 d  e3 E
13. for (i = 0; i < 256; i++)
    
14. {
    
15. input[i] = 0;
    3 |. s6 p7 o( m& {9 r+ [
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
    % A8 h# l1 k% w) a4 W1 R
20.    output：输出结果，高16位是虚部，低16位是实部。
    
21.    input ：输入数据，高16位是虚部，低16位是实部。
    - n8 o; h& K, U. e; x3 l
22.    第三个参数必须是1024。
    
23. */
    
24. cr4_fft_256_stm32(output, input, 256);
    
25. /* 求幅值 */
    
26. PowerMag(256);
    
27. /* 打印输出结果 */
    1 D; L0 ?1 @/ H+ A
28. for (i = 0; i < 256; i++)
    9 Z7 W! P3 `1 E8 t
29. {
    
30. printf("%d\r\n", Mag[i]);
    * U5 W: s7 l, l
31. }
    ' p# v, T: k6 Q7 l3 t, o( `
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+ `7 Q. Q/ |" F5 ]. }& 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计算正弦波中心频率更好些，左右两边的频谱泄露较少，但由于计算的点数较少，Matlb和函数cr4_fft_1024_stm32的效果都比较一般。

9 X" G- E% e* n' {/ H

28.4 函数cr4_fft_64_stm32的使用

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

1. /*
   
2. *********************************************************************************************************
   
3. *        函 数 名: DSP_FFT64
   
4. *        功能说明: 64点FFT实现
   
5. *        形    参：无
   ) Y) S0 T( G+ C9 K
6. *        返 回 值: 无
   $ D1 C! b# y0 S! C5 d& N
7. *********************************************************************************************************
   
8. */
   & x1 X# O) n3 b0 ?1 k/ Z3 w
9. void DSP_FFT64(void)
   : i4 v5 I9 ^( T" H
10. {
    
11. uint16_t i;
    ! {) c  k: x. f' b* r+ Y
12. /* 获得64个采样点 */
    
13. for (i = 0; i < 64; i++)
    
14. {
    
15. input[i] = 0;
    0 Z0 s4 R& ^# [/ u3 s' h4 }  x  P
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. }
    
19. /* 计算64点FFT
    
20.    output：输出结果，高16位是虚部，低16位是实部。
    % B) Y$ P! W. B% a4 j
21.    input ：输入数据，高16位是虚部，低16位是实部。
    
22.    第三个参数必须是1024。
    " Z, F% w" B7 f- W$ P5 k, y" T
23. */
    : s; Y0 T1 |; |0 [
24. cr4_fft_64_stm32(output, input, 64);
    
25. /* 求幅值 */
    + y" f" N- H" S4 c% t; ?
26. PowerMag(64);
    
27. /* 打印输出结果 */
    0 t; m. w8 j* z$ l: l2 b1 K: p
28. for (i = 0; i < 64; i++)
    
29. {
    
30. printf("%d\r\n", Mag[i]);
    
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;          % 真实的频率

3 x4 L; L- B0 m8 L9 A+ X

%波形是由直流分量，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('幅度');

; ]1 y2 B. z/ a: P% c

subplot(2,1,2);

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

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

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

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

28.5 汇编FFT的相频响应求解

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

列出来作为大家的参考。

    相频响应代码求解如下：

1. /*
     Y5 t+ g, W; b4 |4 P" j* w
2. *********************************************************************************************************
   
3. *        函 数 名: Power_Phase_Radians
   $ {! n' s$ y9 a8 N) h
4. *        功能说明: 求相位
   2 T7 M, g* D1 q. I" Z$ I* @
5. *        形    参：_usFFTPoints  FFT点数
   8 q# A% Y9 d3 e
6. *        返 回 值: 无
   
7. *********************************************************************************************************
   
8. */
   $ n4 B- ]2 V: [* F* G& A/ z
9. void Power_Phase_Radians(uint16_t _usFFTPoints)
   
10. {
    9 O  H' C9 s  \2 j) W
11. int16_t lX, lY;
    
12. uint16_t i;
    
13. float32_t phase;
    + Q( t) U; V# L( c
14. for (i=0; i <_usFFTPoints; i++)
    
15. {
    ; H5 Z' |) e+ |2 `& B3 z
16. lX= (output[i]<<16)>>16;  /* 实部 */
    1 n% `! h4 }8 G
17. lY= (output[i] >> 16);    /* 虚部 */
    1 m, v) F" M" x3 S
18.           phase = atan2(lY, lX);    /* atan2求解的结果范围是(-pi, pi], 弧度制 */
    ! J: u( w8 E+ @$ p/ L9 _1 j
19.    
    
20. Phase[i] = phase* 180.0f/3.1415926f;  /* 将求解的结果由弧度转换为角度 */
    ; a) m. [" {, W, x; K* \2 S$ ^1 b
21. }
    
22. }
    
23. 
    8 D9 V- z! O: A
24. /*
    
25. *********************************************************************************************************
    
26. *        函 数 名: DSP_FFTPhase
    
27. *        功能说明: 1024点FFT的相位求解
    
28. *        形    参：无
    
29. *        返 回 值: 无
    9 c$ u- t( g) w" t/ O
30. *********************************************************************************************************
    & Y8 b* E2 r0 q! y
31. */
    / T2 V' [& w6 O" M% E9 @
32. void DSP_FFTPhase(void)
    9 b( L5 J' F6 H9 k; m" O5 E3 w9 S
33. {
    
34. uint16_t i;
    : k6 i$ K2 h9 S' O# i
35. /* 获得1024个采样点 */
    * [5 o0 `/ x( x- W& h% g* {, x0 ^
36. for (i = 0; i < 1024; i++)
    
37. {
    9 M5 U8 ~0 z# O' h/ w" I+ H
38. input[i] = 0;
    
39. /* 波形是由直流分量，50Hz正弦波组成，波形采样率1KHz */
    4 t/ Q" c$ x' O, s" X' G
40. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000 + 3.1415926f/3);
    . O) n  I, E" M& m
41. }
    
42. /* 计算1024点FFT
    8 @1 A4 N' U% ?6 J1 Q8 X+ w' u- E
43.    output：输出结果，高16位是虚部，低16位是实部。
    ! D. p7 }, ]+ s
44.    input ：输入数据，高16位是虚部，低16位是实部。
    
45.    第三个参数必须是1024。
    
46. */
    , ?" H+ u/ C6 |5 z7 U
47. cr4_fft_1024_stm32(output, input, 1024);
    
48. /* 求相频 */
    
49. Power_Phase_Radians(1024);
    
50. /* 打印输出结果 */
    4 X( i) D, L; z2 D
51. for (i = 0; i < 1024; i++)
    * o4 T- P  b, k+ J
52. {
    0 Z. g; l8 I$ V
53. printf("%f\r\n", Phase[i]);
    ) w- c' V% h! ^4 f
54. }
    ; |, t: G2 @0 Q- b+ ^9 @/ I2 F' ?
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;          %真实的频率

$ Q. s9 }7 i" s- D1 |  _- |

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

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

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

) Z3 F! b. M2 f* M$ ?. A" W" ?
' c$ h, X% y5 I9 R7 J$ F; G4 A0 P

subplot(2,1,1);

realvalue = real(y);

imagvalue = imag(y);

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

title('Matlab计算结果');

xlabel('频率');

ylabel('幅度');

8 t/ l/ t. A7 e( G6 M  L( q1 A

subplot(2,1,2);

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

title('STM32计算结果');

xlabel('频率');

ylabel('幅度');

& B! @: S! R1 c- }5 ?5 k

Matlab运行结果如下：

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

2 O5 p6 h4 o( ?  u$ t

28.6总结
4 J$ ~5 ^' u( N本章节主要讲解了汇编FFT的1024点，256点和64点使用方法，有兴趣的可以深入了解汇编代码的实现。

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

标记学习

学习一下。。。

学习一下！！

4 g" D* c* d- ]/ L# ^1 A学习一下。。。收藏了

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