【安富莱DSP教程】第28章 ST官方汇编FFT库应用

席萌0209

特别说明：完整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库说明

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库的移植就完成了。

席萌0209

28.2 函数cr4_fft_1024_stm32的使用

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

1. uint32_t input[1024], output[1024], Mag[1024]; /* 输入，输出和幅值 */
2. /*
3. *********************************************************************************************************
4. *    函 数 名: PowerMag
5. *    功能说明: 求模值
6. *    形    参：_usFFTPoints  FFT点数
7. *    返 回 值: 无
8. *********************************************************************************************************
9. */
10. void PowerMag(uint16_t _usFFTPoints)
11. {
12. int16_t lX,lY;
13. uint16_t i;
14. float32_t mag;
16. /* 计算幅值 */
17. for (i=0; i < _usFFTPoints; i++)
18. {
19.       lX= (output[i]<<16)>>16;          /* 实部*/
20. lY= (output[i]>> 16);             /* 虚部 */   
21. mag = __sqrtf(lX*lX+ lY*lY);      /* 求模 */
22. Mag[i]= mag*2;                    /* 求模后乘以2才是实际模值，直流分量不需要乘2 */
23. }
24.      
25. /* 由于上面多乘了2，所以这里直流分量要除以2 */
26. Mag[0] = Mag[0]>>1;
27. }
29. /*
30. *********************************************************************************************************
31. *    函 数 名: DSP_FFT1024
32. *    功能说明: 1024点FFT实现
33. *    形    参：无
34. *    返 回 值: 无
35. *********************************************************************************************************
36. */
37. void DSP_FFT1024(void)
38. {
39. uint16_t i;
40. /* 获得1024个采样点 */
41. for (i = 0; i < 1024; i++)
42. {
43. input[i] = 0;
44. /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率1KHz */
45. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000) + 512*sin(2*3.1415926f*20*i/1000) ;
46. }
47. /* 计算1024点FFT
48.    output：输出结果，高16位是虚部，低16位是实部。
49.    input ：输入数据，高16位是虚部，低16位是实部。
50.    第三个参数必须是1024。
51. */
52. cr4_fft_1024_stm32(output, input, 1024);
53. /* 求幅值 */
54. PowerMag(1024);
55. /* 打印输出结果 */
56. for (i = 0; i < 1024; i++)
57. {
58. printf("%d\\r\\n", Mag[i]);
59. }
60. }

复制代码

运行函数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;          %真实的频率

%波形是由直流分量，50Hz正弦波和20Hz正弦波组成
x = 1024 + 1024*sin(2*pi*50*t) + 512*sin(2*pi*20*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_1024_stm32计算的结果基本是一直的，但是函数cr4_fft_1024_stm32计算的结果频谱泄露稍严重些（具体情况要具体分析）。

席萌0209

28.3 函数cr4_fft_256_stm32的使用


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

1. /*
2. *********************************************************************************************************
3. *    函 数 名: DSP_FFT256
4. *    功能说明: 256点FFT实现
5. *    形    参：无
6. *    返 回 值: 无
7. *********************************************************************************************************
8. */
9. void DSP_FFT256(void)
10. {
11. uint16_t i;
12. /* 获得256个采样点 */
13. for (i = 0; i < 256; i++)
14. {
15. input[i] = 0;
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
20.    output：输出结果，高16位是虚部，低16位是实部。
21.    input ：输入数据，高16位是虚部，低16位是实部。
22.    第三个参数必须是1024。
23. */
24. cr4_fft_256_stm32(output, input, 256);
25. /* 求幅值 */
26. PowerMag(256);
27. /* 打印输出结果 */
28. for (i = 0; i < 256; i++)
29. {
30. printf("%d\\r\\n", Mag[i]);
31. }
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变换

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

席萌0209

28.4 函数cr4_fft_64_stm32的使用


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

1. /*
2. *********************************************************************************************************
3. *    函 数 名: DSP_FFT64
4. *    功能说明: 64点FFT实现
5. *    形    参：无
6. *    返 回 值: 无
7. *********************************************************************************************************
8. */
9. void DSP_FFT64(void)
10. {
11. uint16_t i;
12. /* 获得64个采样点 */
13. for (i = 0; i < 64; i++)
14. {
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. }
19. /* 计算64点FFT
20.    output：输出结果，高16位是虚部，低16位是实部。
21.    input ：输入数据，高16位是虚部，低16位是实部。
22.    第三个参数必须是1024。
23. */
24. cr4_fft_64_stm32(output, input, 64);
25. /* 求幅值 */
26. PowerMag(64);
27. /* 打印输出结果 */
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;          % 真实的频率

%波形是由直流分量，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计算的结果基本是一直的，但是计算的效果都比较差，主要是因为采样点数太少。

席萌0209

28.5 汇编FFT的相频响应求解


    特别说明：这里相频响应求解暂时有问题（与Matlab求解的结果不一致），但是也把这个求解过程
列出来作为大家的参考。
    相频响应代码求解如下：

1. /*
2. *********************************************************************************************************
3. *    函 数 名: Power_Phase_Radians
4. *    功能说明: 求相位
5. *    形    参：_usFFTPoints  FFT点数
6. *    返 回 值: 无
7. *********************************************************************************************************
8. */
9. void Power_Phase_Radians(uint16_t _usFFTPoints)
10. {
11. int16_t lX, lY;
12. uint16_t i;
13. float32_t phase;
14. for (i=0; i <_usFFTPoints; i++)
15. {
16. lX= (output[i]<<16)>>16;  /* 实部 */
17. lY= (output[i] >> 16);    /* 虚部 */
18.       phase = atan2(lY, lX);    /* atan2求解的结果范围是(-pi, pi], 弧度制 */
19.    
20. Phase[i] = phase* 180.0f/3.1415926f;  /* 将求解的结果由弧度转换为角度 */
21. }
22. }
24. /*
25. *********************************************************************************************************
26. *    函 数 名: DSP_FFTPhase
27. *    功能说明: 1024点FFT的相位求解
28. *    形    参：无
29. *    返 回 值: 无
30. *********************************************************************************************************
31. */
32. void DSP_FFTPhase(void)
33. {
34. uint16_t i;
35. /* 获得1024个采样点 */
36. for (i = 0; i < 1024; i++)
37. {
38. input[i] = 0;
39. /* 波形是由直流分量，50Hz正弦波组成，波形采样率1KHz */
40. input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/1000 + 3.1415926f/3);
41. }
42. /* 计算1024点FFT
43.    output：输出结果，高16位是虚部，低16位是实部。
44.    input ：输入数据，高16位是虚部，低16位是实部。
45.    第三个参数必须是1024。
46. */
47. cr4_fft_1024_stm32(output, input, 1024);
48. /* 求相频 */
49. Power_Phase_Radians(1024);
50. /* 打印输出结果 */
51. for (i = 0; i < 1024; i++)
52. {
53. printf("%f\\r\\n", Phase[i]);
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变换

subplot(2,1,1);
realvalue = real(y);
imagvalue = imag(y);
plot(f, atan2(imagvalue, realvalue)*180/pi);
title('Matlab计算结果');
xlabel('频率');
ylabel('幅度');

subplot(2,1,2);
plot(f, sampledata);   %绘制STM32计算的幅频相应
title('STM32计算结果');
xlabel('频率');
ylabel('幅度');

    Matlab运行结果如下：

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

席萌0209

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