Using STM32 DSP library for FFT transformation

*************************************************** *************************************************** *****
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. */
/ * 8
... ​ ​ ​ ​dsp library in mdk project. * Use trigonometric functions to generate sampling points for FFT calculation * When performing FFT test, call the functions in the following order: * dsp_asm_init(); * dsp_asm_test(); */ #include "stm32f10x.h" #include "dsp_asm.h" #include "stm32_dsp.h" #include "table_fft.h" #include #include /* ************************************************** ************************************************** ***** * LOCAL CONSTANTS ****************************************** ************************************************** ************* */ #define PI2 6.28318530717959 // Comment the lines that you don't want to use. // To simulate FFT, comment out the other predefined // here You can also comment them all out and add NPT_XXX items in MDK's project properties->"C/C++"->"Preprocessor Symbols"-"Define:" // However, the disadvantage of this approach is that every time you modify the XXX data, MDK will compile all files the next time it compiles, which is too slow. //#define NPT_64 64 #define NPT_256 256 //#define NPT_1024 1024 // N=64,Fs/N=50Hz,Max(Valid)=1600Hz // 64-point FFt, sampling rate 3200Hz, frequency resolution 50Hz, measurement Maximum valid frequency 1600Hz #ifdef NPT_64 #define NPT 64 #define Fs 3200 #endif // N=256,Fs/N=25Hz,Max(Valid)=3200Hz // 256-point FFt, sampling rate 6400Hz, frequency resolution 25Hz, The maximum effective frequency of the measurement is 3200Hz #ifdef NPT_256 #define NPT 256 #define Fs 6400 #endif // N=1024,Fs/N=5Hz,Max(Valid)=2560Hz
















































// 1024 point FFt, sampling rate 5120Hz, frequency resolution 5Hz, maximum effective frequency measurement 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
* Simulate sampling data. The sampling data contains 3 frequency sine waves: 50Hz, 2500Hz, 2550Hz
* In the lBUFIN array, the real part of the sampling data is stored in the high word (high 16 bits) of each unit data, and the imaginary part of the sampling data (always 0) is stored in the low word (low 16 bits)
**************************************************************************************************
*/
void dsp_asm_init()
{
  u16 i=0;
  float fx;
  for(i=0;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
* Perform FFT transformation and calculate the amplitude of each harmonic
******************************************************************************************************
*/
void dsp_asm_test()
{
// Select the appropriate FFT function according to the predefined
#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

  // Calculate the amplitude
  dsp_asm_powerMag();
 
// printf("No. Freq Power\n");
// for(i=0;i // {
// printf("%4d,%4d,%10d,%1 0d,%10d\n",i,(u16)((float)i*Fs/NPT),lBUFMAG[i],(lBUFOUT[i]>>16),(lBUFOUT[i]&0xffff));
// }
// printf("*********END**********\r\n");
}
/*
**************************************************************************************************
* Calculate powermag
* Calculate the amplitude of each harmonic
* First decompose lBUFOUT into the real part (X) and the imaginary part (Y), then calculate the assignment (sqrt(X*X+Y*Y)
*********************************************************************************************************
*/
void dsp_asm_powerMag(void)
{
  s16 lX,lY;
  u32 i;
  for(i=0;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);
    }
  }
}


// The author uses the Taurus development board, the CPU is STM32F107VC; JLink V8, MDK-ARM 4.10

// Note the symmetry of the FFT calculation results, that is, for the 256-point calculation results, only the first 128 points of data are valid.
// 64-point FFT calculation result diagram (partial):

 In the figure above, the harmonic frequency corresponding to the array subscript X is: N×Fs/64=N×3200/64=N*50Hz.

lBUFMAG[1] corresponds to the 50Hz harmonic amplitude

In the above figure, due to the FFT resolution of 50HZ, only 1600Hz harmonics can be recognized at most, resulting in incorrect data in the result.
// 256-point FFT calculation result diagram (partial):

 In the figure above, the harmonic frequency corresponding to the array subscript X is: N×Fs/256=N×6400/256=N*25Hz.

lBUFMAG[2] corresponds to 2×25 =50Hz harmonic amplitude

lBUFMAG[100] corresponds to 100×25=2500Hz harmonic amplitude

lBUFMAG[102] corresponds to 102×25=2550Hz harmonic amplitude

// 1024-point FFT calculation result diagram (partial):

 In the figure above, the harmonic frequency corresponding to the array subscript X is: N×Fs/1024=N×5120/1024=N*5Hz.

lBUFMAG[10] corresponds to 10×5 =50Hz harmonic amplitude

lBUFMAG[500] corresponds to 500×5=2500Hz harmonic amplitude

lBUFMAG[510] corresponds to 510×5=2550Hz harmonic amplitude

The analog signal source in this project is: 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs)

The signal is a mixed signal of 1 50Hz, 1 2500Hz, and 1 2550Hz sine waves, all with an amplitude of 4000.