使用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:[e-mail]
* 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
#include


/*
*********************************************************************************************************
* LOCAL CONSTANTS
*********************************************************************************************************
*/
#define PI2 6.28318530717959
// Comment the lines that you dont 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 {
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// {
// 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 {
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);
}
}
}


// 筆者使用的是金牛開發板,CPU為STM32F107VC；JLink V8，MDK-ARM 4.10

本文引用地址：http://www.eepw.com.cn/article/201611/315882.htm

// 注意FFT運算結果的對稱性，也即256點的運算結果，只有前面128點的數據是有效可用的。
// 64點FFT運算結果圖（局部）：

上圖中，數組下標X對應的諧波頻率為：N×Fs/64=N×3200/64=N*50Hz.

lBUFMAG[1] 對應 50Hz諧波幅值

上圖中由於FFT解析度50HZ，最大只能識別1600Hz諧波，導致結果中出現錯誤的數據。
// 256點FFT運算結果圖（局部）：

上圖中，數組下標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運算結果圖（局部）：

上圖中，數組下標X對應的諧波頻率為：N×Fs/1024=N×5120/1024=N*5Hz.

lBUFMAG[10] 對應 10×5 =50Hz諧波幅值

lBUFMAG[500] 對應 500×5=2500Hz諧波幅值

lBUFMAG[510] 對應 510×5=2550Hz諧波幅值

該工程中模擬信號源為：4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs)

信號為1個50Hz、1個2500Hz、1個2550Hz的正弦波混合信號，幅值為均為4000。