/*
*********************************************************************************************************
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 = ((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 ");
//  for(i=0;i<NPT/2;i++)
//  {
//    printf("%4d,%4d,%10d,%10d,%10d ",i,(u16)((float)i*Fs/NPT),lBUFMAG,(lBUFOUT>>16),(lBUFOUT&0xffff));
//  }
//  printf("*********END********** ");
}
/*
*********************************************************************************************************
*   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 << 16) >> 16;
    lY  = (lBUFOUT >> 16);
    {
    float X    = NPT * ((float)lX) /32768;
    float Y    = NPT * ((float)lY) /32768;
    float Mag = sqrt(X*X + Y*Y)/NPT;
    lBUFMAG    = (u32)(Mag * 65536);
    }
  }
}


// 笔者使用的是金牛开发板,CPU为STM32F107VC；JLink V8，MDK-ARM 4.10

// 注意FFT运算结果的对称性，也即256点的运算结果，只有前面128点的数据是有效可用的。
// 64点FFT运算结果图（局部）：

FFT_64 (原文件名:1.JPG)

上图中，数组下标X对应的谐波频率为：N×Fs/64=N×3200/64=N*50Hz.

lBUFMAG[1] 对应 50Hz谐波幅值

上图中由于FFT分辨率50HZ，最大只能识别1600Hz谐波，导致结果中出现错误的数据。
// 256点FFT运算结果图（局部）：


FFT_256 (原文件名:2.JPG)


上图中，数组下标X对应的谐波频率为：N×Fs/256=N×6400/256=N*25Hz.

lBUFMAG[2] 对应 2×25 =50Hz谐波幅值

lBUFMAG[100] 对应 100×25=2500Hz谐波幅值

lBUFMAG[102] 对应 102×25=2550Hz谐波幅值


// 1024点FFT运算结果图（局部）：


FFT_1024 (原文件名:3.JPG)

上图中，数组下标X对应的谐波频率为：N×Fs/1024=N×5120/1024=N*5Hz.

lBUFMAG[10] 对应 10×5 =50Hz谐波幅值

lBUFMAG[500] 对应 500×5=2500Hz谐波幅值

lBUFMAG[510] 对应 510×5=2550Hz谐波幅值



FFT_Code示例代码FFT2.zip(文件大小:490K) (原文件名:FFT2.zip)

更正：
文件dsp_g2.c中dsp_asm_powerMag函数应按如下方式更改：
void dsp_asm_powerMag(void)
{
  s16 lX,lY;
  u32 i;
  for(i=0;i<NPT/2;i++)
  {
    lX  = (lBUFOUT << 16) >> 16;
    lY  = (lBUFOUT >> 16);
        {
        float X        = NPT * ((float)lX) /32768;
        float Y        = NPT * ((float)lY) /32768;
  float Mag = sqrt(X*X + Y*Y)/NPT;
  if(i==0)
    lBUFMAG        = (u32)(Mag * 32768);
  else
          lBUFMAG        = (u32)(Mag * 65536);
    }
  }
}
原因在于计算幅值时直流分量的幅值为sqrt(X*X+Y*Y)/NPT,不用除2