STM32F407的DSP教程第29章 STM32F407移植汇编定点FFT库（64点，256点和1024点）

Posted 2021-06-10 安富莱电子

完整版教程下载地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547

第29章       STM32F407移植汇编定点FFT库（64点，256点和1024点）

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

29.1 汇编FFT库说明

29.2 函数cr4_fft_1024_stm32的使用(含幅频和相频响应)

29.3 函数cr4_fft_256_stm32的使用

29.4 函数cr4_fft_64_stm32的使用

29.5 实验例程说明(MDK)

29.6 实验例程说明(IAR)

29.7 总结。

 

 

29.1 汇编FFT库说明

29.1.1 描述

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

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

•   N得满足4ⁿ（n =1，2, 3…..），也就是以4为基数。
•   采样信号必须是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。

29.1.2 汇编库的移植

注：这里以MDK为例进行说明，IAR是一样的。

这个汇编库的移植比较简单，从本章配套例子User文件夹复制fft文件夹到自己的工程：

 

注意路径\\User\\fft\\src\\asm下有三个文件夹，分布是arm，gcc和iar，其中arm可用于MDK，gcc可用于Embedded Studio，iar可用于IAR FOR ARM。

 

三个文件夹里面都是如下几个文件，只是用于不用的编译器：

 

然后把FFT源文件的三个FFT汇编文件和两个头文件添加上即可，添加后效果如下(注意不同编译器添相应汇编文件)：

 

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

 

最后别忘了添加路径：

 

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

29.2 函数cr4_fft_1024_stm32的使用(含幅频和相频响应)

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

uint32_t input[1024], output[1024], Mag[1024];/* 输入，输出和幅值 */ 
float32_t Phase[1024]; /* 相位*/

/*
*********************************************************************************************************
*    函 数 名: PowerMag
*    功能说明: 求模值
*    形    参：_usFFTPoints  FFT点数
*    返 回 值: 无
*********************************************************************************************************
*/
void PowerMag(uint16_t _usFFTPoints)        
{
    int16_t lX,lY;
    uint16_t i;
    float32_t mag;

    /* 计算幅值 */
    for (i=0; i < _usFFTPoints; i++)
    {
         lX= (output[i]<<16)>>16;          /* 实部*/
        lY= (output[i]>> 16);             /* 虚部 */    
        arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag);   /* 求模 */
        Mag[i]= mag*2;                    /* 求模后乘以2才是实际模值，直流分量不需要乘2 */
    }
     
    /* 由于上面多乘了2，所以这里直流分量要除以2 */
    Mag[0] = Mag[0]>>1;
}

/*
*********************************************************************************************************
*    函 数 名: Power_Phase_Radians
*    功能说明: 求相位
*    形    参：_usFFTPoints  FFT点数， uiCmpValue 阀值
*    返 回 值: 无
*********************************************************************************************************
*/
void Power_Phase_Radians(uint16_t _usFFTPoints, uint32_t _uiCmpValue)        
{
    int16_t lX, lY;
    uint16_t i;
    float32_t phase;
    float32_t mag;
    
    
    for (i=0; i <_usFFTPoints; i++)
    {
        lX= (output[i]<<16)>>16;  /* 实部 */
        lY= (output[i] >> 16);    /* 虚部 */ 
        
         phase = atan2(lY, lX);            /* atan2求解的结果范围是(-pi, pi], 弧度制 */
        arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag);   /* 求模 */
        
        if(_uiCmpValue > mag)
        {
            Phase[i] = 0;            
        }
        else
        {
            Phase[i] = phase* 180.0f/3.1415926f;  /* 将求解的结果由弧度转换为角度 */
        }
    }
}

/*
*********************************************************************************************************
*    函 数 名: DSP_FFTPhase
*    功能说明: 1024点FFT的相位求解
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
void DSP_FFTPhase(void)
{
    uint16_t i;
    
    /* 获得1024个采样点 */
    for (i = 0; i < 1024; i++)
    {
        /* 波形是由直流分量，50Hz正弦波组成，波形采样率1024 */
        input[i] = 1024 + 1024*cos(2*3.1415926f*50*i/1024 + 3.1415926f/3);
    }
    
    /* 计算1024点FFT 
       output：输出结果，高16位是虚部，低16位是实部。
       input ：输入数据，高16位是虚部，低16位是实部。
       第三个参数必须是1024。
    */
    cr4_fft_1024_stm32(output, input, 1024); 
    
    /* 求幅值 */
    PowerMag(1024);
    
    /* 打印输出结果 */
    for (i = 0; i < 1024; i++)
    {
        printf("%d\\r\\n", Mag[i]);
    }
    
    printf("=========================================\\r\\n");    
    
    /* 求相频 */
    Power_Phase_Radians(1024, 100);
    
    
    /* 打印输出结果 */
    for (i = 0; i < 1024; i++)
    {
        printf("%f\\r\\n", Phase[i]);
    }
}

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

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

Fs = 1024;               % 采样率
N  = 1024;               % 采样点数
n  = 0:N-1;              % 采样序列
t  = 0:1/Fs:1-1/Fs;      % 时间序列
f = n * Fs / N;          %真实的频率

%波形是由直流分量，50Hz正弦波正弦波组成
x = 1024 + 1024*cos(2*pi*50*t + pi/3)   ;  
y = fft(x, N);               %对原始信号做FFT变换
Mag = abs(y);

subplot(2,2,1);
MagAct = Mag *2  / N;
MagAct(1) = MagAct(1)/2;
plot(f, MagAct); 
title(\'Matlab计算幅频响应\');
xlabel(\'频率\');
ylabel(\'赋值\');

subplot(2,2,2);
realvalue = real(y);
imagvalue = imag(y);
plot(f, atan2(imagvalue, realvalue)*180/pi.*(Mag>=1024*20)); 
title(\'Matlab计算相频响应\');
xlabel(\'频率\');
ylabel(\'相角\');

subplot(2,2,3);
plot(f, sampledata1);  %绘制STM32计算的幅频相应
title(\'STM32计算幅频响应\');
xlabel(\'频率\');
ylabel(\'赋值\');

subplot(2,2,4);
plot(f, sampledata2);   %绘制STM32计算的相频相应
title(\'STM32计算相频响应\');
xlabel(\'频率\');
ylabel(\'相角\');

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

 

从上面的对比结果中可以看出，Matlab和函数cr4_fft_1024_stm32计算的结果基本是一直的。幅频响应求出的幅值和相频响应中的求出的初始相角都是没问题的。

29.3 函数cr4_fft_256_stm32的使用

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

/*
*********************************************************************************************************
*    函 数 名: DSP_FFT256
*    功能说明: 256点FFT实现
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
void DSP_FFT256(void)
{
    uint16_t i;
    
    /* 获得256个采样点 */
    for (i = 0; i < 256; i++)
    {
        /* 波形是由直流分量，50Hz正弦波和20Hz正弦波组成，波形采样率200Hz */
        input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/200) + 512*sin(2*3.1415926f*20*i/200) ;
    }
    
    /* 计算256点FFT 
       output：输出结果，高16位是虚部，低16位是实部。
       input ：输入数据，高16位是虚部，低16位是实部。
       第三个参数必须是1024。
    */
    cr4_fft_256_stm32(output, input, 256); 
    
    /* 求幅值 */
    PowerMag(256);
    
    /* 打印输出结果 */
    for (i = 0; i < 256; i++)
    {
        printf("%d\\r\\n", Mag[i]);
    }
}

运行函数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);
MagAct = Mag *2 / N;
MagAct(1) = MagAct(1)/2;
plot(f, MagAct);          %绘制幅频相应曲线
title(\'Matlab计算结果\');
xlabel(\'频率\');
ylabel(\'幅度\');

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

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

 

从上面的对比结果中可以看出，Matlab和函数cr4_fft_256_stm32计算的结果基本是一直的，但频率泄露略多。

29.4 函数cr4_fft_64_stm32的使用

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

/*
*********************************************************************************************************
*    函 数 名: DSP_FFT64
*    功能说明: 64点FFT实现
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
void DSP_FFT64(void)
{
    uint16_t i;
    
    /* 获得64个采样点 */
    for (i = 0; i < 64; i++)
    {
        /* 波形是由直流分量，5Hz正弦波和10Hz正弦波组成，波形采样率60Hz */
        input[i] = 1024 + 1024*sin(2*3.1415926f*5*i/60) + 512*sin(2*3.1415926f*10*i/60) ;
    }
    
    /* 计算64点FFT 
       output：输出结果，高16位是虚部，低16位是实部。
       input ：输入数据，高16位是虚部，低16位是实部。
       第三个参数必须是1024。
    */
    cr4_fft_64_stm32(output, input, 64); 
    
    /* 求幅值 */
    PowerMag(64);
    
    /* 打印输出结果 */
    for (i = 0; i < 64; i++)
    {
        printf("%d\\r\\n", Mag[i]);
    }
}

运行函数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);
MagAct = Mag *2 / N;
MagAct(1) = MagAct(1)/2;
plot(f, MagAct);          %绘制幅频相应曲线
title(\'Matlab计算结果\');
xlabel(\'频率\');
ylabel(\'幅度\');

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

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

 

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

29.5 实验例程说明（MDK）

配套例子：

V5-219_STM32F407移植ST汇编定点FFT库(64点，256点和1024点)

实验目的：

1. 学习ST汇编定点FFT库（64点，256点和1024点）

实验内容：

1. 启动一个自动重装软件定时器，每100ms翻转一次LED2。
2. 按下按键K1，串口打印1024点FFT的幅频响应和相频响应。
3. 按下按键K2，串口打印256点FFT的幅频响应。
4. 按下按键K3，串口打印64点FFT的幅频响应。

使用AC6注意事项

特别注意附件章节C的问题

上电后串口打印的信息：

波特率 115200，数据位 8，奇偶校验位无，停止位 1。

 

RTT方式打印信息：

 

程序设计：

  系统栈大小分配：

 

  硬件外设初始化

硬件外设的初始化是在 bsp.c 文件实现：

/*
*********************************************************************************************************
*    函 数 名: bsp_Init
*    功能说明: 初始化所有的硬件设备。该函数配置CPU寄存器和外设的寄存器并初始化一些全局变量。只需要调用一次
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
void bsp_Init(void)
{
    /* 
       STM32F407 HAL 库初始化，此时系统用的还是F407自带的16MHz，HSI时钟:
       - 调用函数HAL_InitTick，初始化滴答时钟中断1ms。
       - 设置NVIV优先级分组为4。
     */
    HAL_Init();

    /* 
       配置系统时钟到168MHz
       - 切换使用HSE。
       - 此函数会更新全局变量SystemCoreClock，并重新配置HAL_InitTick。
    */
    SystemClock_Config();

    /* 
       Event Recorder：
       - 可用于代码执行时间测量，MDK5.25及其以上版本才支持，IAR不支持。
       - 默认不开启，如果要使能此选项，务必看V5开发板用户手册第8章
    */    
#if Enable_EventRecorder == 1  
    /* 初始化EventRecorder并开启 */
    EventRecorderInitialize(EventRecordAll, 1U);
    EventRecorderStart();
#endif
    
    bsp_InitKey();        /* 按键初始化，要放在滴答定时器之前，因为按钮检测是通过滴答定时器扫描 */
    bsp_InitTimer();      /* 初始化滴答定时器 */
    bsp_InitUart();    /* 初始化串口 */
    bsp_InitExtIO();   /* 初始化扩展IO */
    bsp_InitLed();        /* 初始化LED */        
}

  主功能：

主程序实现如下操作：

•   启动一个自动重装软件定时器，每100ms翻转一次LED2。
•   按下按键K1，串口打印1024点FFT的幅频响应和相频响应。
•   按下按键K2，串口打印256点FFT的幅频响应。
•   按下按键K3，串口打印64点FFT的幅频响应。

/*
*********************************************************************************************************
*    函 数 名: main
*    功能说明: c程序入口
*    形    参: 无
*    返 回 值: 错误代码(无需处理)
*********************************************************************************************************
*/
int main(void)
{
    uint8_t ucKeyCode;        /* 按键代码 */
    

    bsp_Init();        /* 硬件初始化 */
    PrintfLogo();    /* 打印例程信息到串口1 */

    PrintfHelp();    /* 打印操作提示信息 */
    

    bsp_StartAutoTimer(0, 100);    /* 启动1个100ms的自动重装的定时器 */

    /* 进入主程序循环体 */
    while (1)
    {
        bsp_Idle();        /* 这个函数在bsp.c文件。用户可以修改这个函数实现CPU休眠和喂狗 */

        if (bsp_CheckTimer(0))    /* 判断定时器超时时间 */
        {
            /* 每隔500ms 进来一次 */
            bsp_LedToggle(4);    /* 翻转LED2的状态 */   
        }
        
        ucKeyCode = bsp_GetKey();    /* 读取键值, 无键按下时返回 KEY_NONE = 0 */
        if (ucKeyCode != KEY_NONE)
        {
            switch (ucKeyCode)
            {
                case KEY_DOWN_K1:            /* K1键按下 */
                    DSP_FFTPhase();
                    break;
                
                case KEY_DOWN_K2:            /* K2键按下 */
                    DSP_FFT256();
                    break;
                
                case KEY_DOWN_K3:            /* K3键按下 */
                    DSP_FFT64();
                    break;
                
                    
                default:
                    /* 其它的键值不处理 */
                    break;
            }
        }

    }
}

29.6 实验例程说明（IAR）

配套例子：

V5-219_STM32F407移植ST汇编定点FFT库(64点，256点和1024点)

实验目的：

1. 学习ST汇编定点FFT库（64点，256点和1024点）

实验内容：

1. 启动一个自动重装软件定时器，每100ms翻转一次LED2。
2. 按下按键K1，串口打印1024点FFT的幅频响应和相频响应。
3. 按下按键K2，串口打印256点FFT的幅频响应。
4. 按下按键K3，串口打印64点FFT的幅频响应。

上电后串口打印的信息：

波特率 115200，数据位 8，奇偶校验位无，停止位 1。

 

RTT方式打印信息：

 

程序设计：

  系统栈大小分配：

 

  硬件外设初始化

硬件外设的初始化是在 bsp.c 文件实现：

/*
*********************************************************************************************************
*    函 数 名: bsp_Init
*    功能说明: 初始化所有的硬件设备。该函数配置CPU寄存器和外设的寄存器并初始化一些全局变量。只需要调用一次
*    形    参：无
*    返 回 值: 无
*********************************************************************************************************
*/
void bsp_Init(void)
{
    /* 
       STM32F407 HAL 库初始化，此时系统用的还是F407自带的16MHz，HSI时钟:
       - 调用函数HAL_InitTick，初始化滴答时钟中断1ms。
       - 设置NVIV优先级分组为4。
     */
    HAL_Init();

    /* 
       配置系统时钟到168MHz
       - 切换使用HSE。
       - 此函数会更新全局变量SystemCoreClock，并重新配置HAL_InitTick。
    */
    SystemClock_Config();

    /* 
       Event Recorder：
       - 可用于代码执行时间测量，MDK5.25及其以上版本才支持，IAR不支持。
       - 默认不开启，如果要使能此选项，务必看V5开发板用户手册第8章
    */    
#if Enable_EventRecorder == 1  
    /* 初始化EventRecorder并开启 */
    EventRecorderInitialize(EventRecordAll, 1U);
    EventRecorderStart();
#endif
    
    bsp_InitKey();        /* 按键初始化，要放在滴答定时器之前，因为按钮检测是通过滴答定时器扫描 */
    bsp_InitTimer();      /* 初始化滴答定时器 */
    bsp_InitUart();    /* 初始化串口 */
    bsp_InitExtIO();   /* 初始化扩展IO */
    bsp_InitLed();        /* 初始化LED */        
}

  主功能：

主程序实现如下操作：

•   启动一个自动重装软件定时器，每100ms翻转一次LED2。
•   按下按键K1，串口打印1024点FFT的幅频响应和相频响应。
•   按下按键K2，串口打印256点FFT的幅频响应。
•   按下按键K3，串口打印64点FFT的幅频响应。

/*
*********************************************************************************************************
*    函 数 名: main
*    功能说明: c程序入口
*    形    参: 无
*    返 回 值: 错误代码(无需处理)
*********************************************************************************************************
*/
int main(void)
{
    uint8_t ucKeyCode;        /* 按键代码 */
    

    bsp_Init();        /* 硬件初始化 */
    PrintfLogo();    /* 打印例程信息到串口1 */

    PrintfHelp();    /* 打印操作提示信息 */
    

    bsp_StartAutoTimer(0, 100);    /* 启动1个100ms的自动重装的定时器 */

    /* 进入主程序循环体 */
    while (1)
    {
        bsp_Idle();        /* 这个函数在bsp.c文件。用户可以修改这个函数实现CPU休眠和喂狗 */

        if (bsp_CheckTimer(0))    /* 判断定时器超时时间 */
        {
            /* 每隔500ms 进来一次 */
            bsp_LedToggle(4);    /* 翻转LED2的状态 */   
        }
        
        ucKeyCode = bsp_GetKey();    /* 读取键值, 无键按下时返回 KEY_NONE = 0 */
        if (ucKeyCode != KEY_NONE)
        {
            switch (ucKeyCode)
            {
                case KEY_DOWN_K1:            /* K1键按下 */
                    DSP_FFTPhase();
                    break;
                
                case KEY_DOWN_K2:            /* K2键按下 */
                    DSP_FFT256();
                    break;
                
                case KEY_DOWN_K3:            /* K3键按下 */
                    DSP_FFT64();
                    break;
                
                    
                default:
                    /* 其它的键值不处理 */
                    break;
            }
        }

    }
}

29.7 总结

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