使用STM32 的DSP库进行FFT变换

/*

*********************************************************************************************************

FileName:dsp_asm.h

*********************************************************************************************************

*/



#ifndef  __DSP_ASM_H__

#define  __DSP_ASM_H__

*********************************************************************************************************

*                                            FUNCTION PROTOTYPES

*********************************************************************************************************

*/



void    dsp_asm_test(void);

void  dsp_asm_init(void);



#endif                                                          /* End of module include.                               */

/*8888888888888888888888888888888888888888888888888888888888888888*/

/*8888888888888888888888888888888888888888888888888888888888888888*/

/*

* FileName:dsp_asm.c

* Author:Bobby.Chen

* Email:heroxx@163.com

* Date:2010-08-11

* Description:This file showes how to use the dsp library in mdk project.

*                  使用三角函数生成采样点，供FFT计算

*                  进行FFT测试时，按下面顺序调用函数即可：

*                   dsp_asm_init();

*                   dsp_asm_test();

*/

#include "stm32f10x.h"

#include "dsp_asm.h"

#include "stm32_dsp.h"

#include "table_fft.h"

#include <stdio.h>

#include <math.h>





/*

*********************************************************************************************************

*                                           LOCAL CONSTANTS

*********************************************************************************************************

*/

#define PI2 6.28318530717959

// Comment the lines that you don't want to use.

// 要模拟FFT，请注释掉其他的预定义

// 此处也可以全部注释掉，在MDK的工程属性->"C/C++"->"Preprocessor Symbols"-"Define:"中添加NPT_XXX项目

// 但是这样做法的缺点是每次修改XXX数据，都会导致MDK下次编译时会编译全部文件，速度太慢。

//#define NPT_64 64

#define NPT_256 256

//#define NPT_1024 1024



// N=64,Fs/N=50Hz,Max(Valid)=1600Hz

// 64点FFt，采样率3200Hz,频率分辨率50Hz，测量最大有效频率1600Hz

#ifdef NPT_64

#define NPT 64

#define Fs  3200

#endif



// N=256,Fs/N=25Hz,Max(Valid)=3200Hz

// 256点FFt，采样率6400Hz,频率分辨率25Hz，测量最大有效频率3200Hz

#ifdef NPT_256

#define NPT 256

#define Fs  6400

#endif



// N=1024,Fs/N=5Hz,Max(Valid)=2560Hz

// 1024点FFt，采样率5120Hz,频率分辨率5Hz，测量最大有效频率2560Hz

#ifdef NPT_1024

#define NPT 1024

#define Fs  5120

#endif



/*

*********************************************************************************************************

*                                       LOCAL GLOBAL VARIABLES

*********************************************************************************************************

*/

extern  uint16_t  TableFFT[];

long lBUFIN[NPT];         /* Complex input vector */

long lBUFOUT[NPT];        /* Complex output vector */

long lBUFMAG[NPT];/* Magnitude vector */

/*

*********************************************************************************************************

*                                      LOCAL FUNCTION PROTOTYPES

*********************************************************************************************************

*/

void dsp_asm_powerMag(void);



/*

*********************************************************************************************************

*   Initialize data tables for lBUFIN

* 模拟采样数据,采样数据中包含3种频率正弦波：50Hz，2500Hz，2550Hz

* lBUFIN数组中，每个单元数据高字(高16位)中存储采样数据的实部，低字(低16位)存储采样数据的虚部(总是为0)

*********************************************************************************************************

*/

void  dsp_asm_init()

{

  u16 i=0;

  float fx;

  for(i=0;i<NPT;i++)

  {

    fx  = 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs);

    lBUFIN[i] = ((s16)fx)<<16;

  }

}



/*

*********************************************************************************************************

*   Test FFT,calculate powermag

*   进行FFT变换，并计算各次谐波幅值

*********************************************************************************************************

*/

void  dsp_asm_test()

{

// 根据预定义选择合适的FFT函数

#ifdef NPT_64

  cr4_fft_64_stm32(lBUFOUT, lBUFIN, NPT);

#endif



#ifdef NPT_256

  cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT);

#endif



#ifdef NPT_1024

  cr4_fft_1024_stm32(lBUFOUT, lBUFIN, NPT);

#endif



  // 计算幅值

  dsp_asm_powerMag();

 

//  printf("No. Freq      Power\n");

//  for(i=0;i<NPT/2;i++)

//  {

//    printf("%4d,%4d,%10d,%10d,%10d\n",i,(u16)((float)i*Fs/NPT),lBUFMAG[i],(lBUFOUT[i]>>16),(lBUFOUT[i]&0xffff));

//  }

//  printf("*********END**********\r\n");

}

/*

*********************************************************************************************************

*   Calculate powermag

*   计算各次谐波幅值

*   先将lBUFOUT分解成实部(X)和虚部(Y)，然后计算赋值(sqrt(X*X+Y*Y)

*********************************************************************************************************

*/

void dsp_asm_powerMag(void)

{

  s16 lX,lY;

  u32 i;

  for(i=0;i<NPT/2;i++)

  {

    lX  = (lBUFOUT[i] << 16) >> 16;

    lY  = (lBUFOUT[i] >> 16);

    {

    float X    = NPT * ((float)lX) /32768;

    float Y    = NPT * ((float)lY) /32768;

    float Mag = sqrt(X*X + Y*Y)/NPT;

    lBUFMAG[i]    = (u32)(Mag * 65536);

    }

  }

}