前些天理解了fft變換的理論和對其工程應用進行了實例分析，詳見我的名為《C語言實現fft理論基礎與工程應用的實例分析》的博客，用C語言編寫的fft算法比較容易看懂，但帶來的缺點就是執行效率低，對於要求實時操作（例如電機控制）的反應速度不夠靈敏。本篇內容將簡要分析STM32自帶的DSP庫文件，其用匯編語言編寫，代碼執行效率明顯優於C語言，ST公司封裝好了了庫文件，我們不必看懂其匯編代碼，只要會調用接口函數即可。
1,代碼分析
首先我們需要在一個已經建立好的工程文件里添加如下編譯路徑：
工程需要添加的文件如下圖：
1 /* 2 ********************************************************************************************************* 3 * MICIRUM BOARD SUPPORT PACKAGE 4 * 5 * (c) Copyright 2007; Micrium, Inc.; Weston, FL 6 * 7 * All rights reserved. Protected by international copyright laws. 8 * Knowledge of the source code may NOT be used to develop a similar product. 9 * Please help us continue to provide the Embedded community with the finest 10 * software available. Your honesty is greatly appreciated. 11 ********************************************************************************************************* 12 */ 13 14 /* 15 ********************************************************************************************************* 16 * 17 * BOARD SUPPORT PACKAGE 18 * 19 * ST Microelectronics STM32 20 * with the 21 * STM3210B-EVAL Evaluation Board 22 * 23 * Filename : bsp.c 24 * Version : V1.00 25 * Programmer(s) : Brian Nagel 26 ********************************************************************************************************* 27 */ 28 29 /* 30 ********************************************************************************************************* 31 * INCLUDE FILES 32 ********************************************************************************************************* 33 */ 34 35 #define DSP_ASM 36 #include "stm32f10x.h" 37 #include "dsp_asm.h" 38 #include "stm32_dsp.h" 39 #include "table_fft.h" 40 #include <stdio.h> 41 #include <math.h> 42 43 44 /* 45 ********************************************************************************************************* 46 * LOCAL CONSTANTS 47 ********************************************************************************************************* 48 */ 49 #define PI2 6.28318530717959 50 #define NPT_1024 1024 51 //#define NPT_256 256 52 //#define NPT_1024 1024 53 54 // N=64,Fs/N=50Hz,Max(Valid)=1600Hz 55 #ifdef NPT_64 56 #define NPT 64 57 #define Fs 3200 58 #endif 59 60 // N=256,Fs/N=25Hz,Max(Valid)=3200Hz 61 #ifdef NPT_256 62 #define NPT 256 63 #define Fs 6400 64 #endif 65 66 // N=1024,Fs/N=5Hz,Max(Valid)=2560Hz 67 #ifdef NPT_1024 68 #define NPT 1024 69 #define Fs 5120 70 #endif 71 72 73 /* 74 ********************************************************************************************************* 75 * LOCAL DATA TYPES 76 ********************************************************************************************************* 77 */ 78 79 80 /* 81 ********************************************************************************************************* 82 * LOCAL TABLES 83 ********************************************************************************************************* 84 */ 85 86 87 /* 88 ********************************************************************************************************* 89 * LOCAL GLOBAL VARIABLES 90 ********************************************************************************************************* 91 */ 92 long lBUFIN[NPT]; /* Complex input vector */ 93 long lBUFOUT[NPT]; /* Complex output vector */ 94 long lBUFMAG[NPT];/* Magnitude vector */ 95 /* 96 ********************************************************************************************************* 97 * LOCAL FUNCTION PROTOTYPES 98 ********************************************************************************************************* 99 */ 100 void dsp_asm_powerMag(void); 101 102 /* 103 ********************************************************************************************************* 104 * LOCAL CONFIGURATION ERRORS 105 ********************************************************************************************************* 106 */ 107 108 109 /* 110 ****************************************************************************************************************************** 111 ****************************************************************************************************************************** 112 ** Global Functions 113 ****************************************************************************************************************************** 114 ****************************************************************************************************************************** 115 */ 116 117 void dsp_asm_init() 118 { 119 u16 i=0; 120 float fx; 121 for(i=0;i<NPT;i++) 122 { 123 fx = 4000 * sin(PI2*i*50.0/Fs) + 4000 * sin(PI2*i*2500.0/Fs) + 4000*sin(PI2*i*2550.0/Fs); 124 lBUFIN[i] = ((long)fx)<<16; 125 } 126 } 127 128 void dsp_asm_test() 129 { 130 131 #ifdef NPT_64 132 cr4_fft_64_stm32(lBUFOUT, lBUFIN, NPT); 133 #endif 134 135 #ifdef NPT_256 136 cr4_fft_256_stm32(lBUFOUT, lBUFIN, NPT); 137 #endif 138 139 #ifdef NPT_1024 140 cr4_fft_1024_stm32(lBUFOUT, lBUFIN, NPT); 141 #endif 142 143 // 計算幅值 144 dsp_asm_powerMag(); 145 146 } 147 148 void dsp_asm_powerMag(void) 149 { 150 s16 lX,lY; 151 u32 i; 152 for(i=0;i<NPT/2;i++) 153 { 154 lX = (lBUFOUT[i] << 16) >> 16; 155 lY = (lBUFOUT[i] >> 16); 156 { 157 float X = NPT * ((float)lX) /32768; 158 float Y = NPT * ((float)lY) /32768; 159 float Mag = sqrt(X*X + Y*Y)/NPT; 160 lBUFMAG[i] = (u32)(Mag * 65536); 161 } 162 } 163 }
1 /* 2 ********************************************************************************************************* 3 * MICIRUM BOARD SUPPORT PACKAGE 4 * 5 * (c) Copyright 2007; Micrium, Inc.; Weston, FL 6 * 7 * All rights reserved. Protected by international copyright laws. 8 * Knowledge of the source code may NOT be used to develop a similar product. 9 * Please help us continue to provide the Embedded community with the finest 10 * software available. Your honesty is greatly appreciated. 11 ********************************************************************************************************* 12 */ 13 14 /* 15 ********************************************************************************************************* 16 * 17 * BOARD SUPPORT PACKAGE 18 * 19 * ST Microelectronics STM32 20 * with the 21 * STM3210B-EVAL Evaluation Board 22 * 23 * Filename : bsp.h 24 * Version : V1.00 25 * Programmer(s) : Brian Nagel 26 ********************************************************************************************************* 27 */ 28 29 #ifndef __DSP_ASM_H__ 30 #define __DSP_ASM_H__ 31 32 /* 33 ********************************************************************************************************* 34 * EXTERNS 35 ********************************************************************************************************* 36 */ 37 38 #ifdef DSP_ASM 39 #define DSP_EXT 40 #else 41 #define DSP_EXT extern 42 #endif 43 44 /* 45 ********************************************************************************************************* 46 * INCLUDE FILES 47 ********************************************************************************************************* 48 */ 49 50 51 52 53 /* 54 ********************************************************************************************************* 55 * GLOBAL VARIABLES 56 ********************************************************************************************************* 57 */ 58 59 60 /* 61 ********************************************************************************************************* 62 * MACRO'S 63 ********************************************************************************************************* 64 */ 65 66 67 /* 68 ********************************************************************************************************* 69 * FUNCTION PROTOTYPES 70 ********************************************************************************************************* 71 */ 72 73 void dsp_asm_test(void); 74 void dsp_asm_init(void); 75 76 #endif /* End of module include. */
着重分析下dsp_asm_powerMag()函數的作用，其函數就是求幅值，首先定義的的一個16位的有符號的數據IX 和IY 這兩個只是中間變量，然后定義的i，是32位的無符號型。語句的目的是Mag = sqrt(X*X + Y*Y)/NPT。但直接這么寫不符合DSP的計算習慣也就是不符合浮點運算的習慣。因此語句在for函數i寫道 lX = (lBUFOUT[i] << 16) >> 16 就是取32位的i的低16位數據，lY = (lBUFOUT[i] >> 16);是取高16位數據。下面的兩句
float X = NPT * ((float)lX) /32768;
float Y = NPT * ((float)lY) /32768
目的就是把數據浮點化，至於為什么是除以32768 。可以這么說，浮點化就好像10進制里面的科學計數法。32768=2的15次。除以32768也就是去除了浮點數后面的那個基數，只剩下前面的。比如1991 改寫成1.991*10的三次冪，再除以10的三次方，只剩下1.991，便於余下的運算。至於最后一句要乘以65536是因為我們定義的數據和我們需要求得的數據都是無符號32位的，之前已經把32位的數據拆開又分別浮點化了又開了個根號，所以再把它變回來 只需要乘以2的16次，也就是65536.比如說問你什么時候生日，你說是19911030，然而DSP是不習慣這么干的，他需要把它拆開為1991和1030。再寫成1.991x10的3次方和1.030x10的3次方。然后才能進行其他的運算。
這里是ST公司采用了DSP專用芯片（主要是指TI）的寫法，也就是說盡管DSP的芯片類型很多，數據變量的定義也各有差異，但原理是一樣的，最終還是要采用DSP習慣的運算方式。至於為什么一定要采用浮點運算，因為機器是傻子，然而TI公司的工程師是天才。
main函數中我們只需在while（1）前加上dsp_asm_init(); dsp_asm_test();即可。
2，實驗現象
注意FFT運算結果的對稱性，也即256點的運算結果，只有前面128點的數據是有效可用的。
① N=64,Fs/N=50Hz,Max(Valid)=1600Hz，64點FFt，采樣率3200Hz,頻率分辨率50Hz，測量最大有效頻率1600Hz
64點FFT運算結果圖（局部）：
上圖中，數組下標X對應的諧波頻率為：N×Fs/64=N×3200/64=N*50Hz.
lBUFMAG[1] 對應 50Hz諧波幅值。
上圖中由於FFT分辨率50HZ，最大只能識別1600Hz諧波，導致結果中出現錯誤的數據。
②N=256,Fs/N=25Hz,Max(Valid)=3200Hz，256點FFt，采樣率6400Hz,頻率分辨率25Hz，測量最大有效頻率3200Hz
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諧波幅值
③N=1024,Fs/N=5Hz,Max(Valid)=2560Hz，1024點FFt，采樣率5120Hz,頻率分辨率5Hz，測量最大有效頻率2560Hz
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。
3，工程源碼下載地址
LL STM32-DSP-fft.rar
115網盤禮包碼：5lbdd57l2hkq
http://115.com/lb/5lbdd57l2hkq