STM32使用FFT變換計算THD（20年四川省電子設計大賽E題軟件部分）

注： 本篇內容意在使不理解FFT變換的讀者也可以通過使用FFT來計算總諧波失真

FFT變換

根據總諧波失真的定義：

\[THD = \frac{\sqrt{\sum_{n=0}^{\infty}{G_{n}^{2}}}}{G_0} (G_0為基波，G_n 為高次諧波) \]

可知，要計算THD需要知道基波分量和各個諧波分量的大小。

​ FFT也叫快速傅里葉變換，是離散時間傅里葉變換的一種實現手段，變換的本質還是離散時間傅里葉變換（DTFT）。傅里葉變換可以將信號從時域展開到頻域，通過FFT變換即可實現對信號基波和諧波分量的計算。

STM32F4進行FFT

​ 如果不熟悉傅里葉變換，也無關緊要，你只要知道FFT變換可以得到信號在各個頻率上的分量，以STM32F407VGx為例，我們可以使用ST提供的DSP庫快速實現FFT變換。

​ 這里使用cube MX生成一個簡單的工程作為模板，首先加入一些使用DSP庫所需要的宏定義：

,ARM_MATH_CM4,__CC_ARM,ARM_MATH_MATRIX_CHECK,ARM_MATH_ROUNDING

然后，進入CubeMX的安裝目錄，找到DSP庫下的幾個文件

首先lib文件：

​ arm_cortexM4lf_math.lib

cortexM4 是指 arm的cortexM4內核，后面的l 是指芯片使用小端存放，f是浮點運算，STM32F4 默認是小端存放所以使用這個庫文件

位置C:\Users\Username\STM32Cube\Repository\STM32Cube_FW_F4_V1.25.1\Drivers\CMSIS\Lib\ARM 下、

然后是頭文件

位置在

C:\Users\Username\STM32Cube\Repository\STM32Cube_FW_F4_V1.25.1\Drivers\CMSIS\DSP\Include

將這三個文件加入工程中

注意：

如果是和我一樣通過CubeMX生成的工程，且在code Generator中選擇了Copy all library into the floder 如圖

則可以在工程文件夾下找到所需要的的文件，文件目錄結構與cubeMX下的路徑相同

如果覺得上述操作過於麻煩，

可以使用cubeMX自帶的例程模板，在此模板的基礎上進行改進

DSP庫的Example 位於：安裝目錄\STM32Cube\Repository\STM32Cube_FW_F4_V1.25.1\Drivers\CMSIS\DSP\Examples\

有ARM和GCC兩個版本，這里以ARM為例，

不過此工程並未未添加宏定義，且使用的cortexM0內核，需要做一定的修改。個人不推薦使用此文件構建工程，這個工程更適合與用來熟悉函數的用法和功能。

然后在已經配置好的工程中編寫代碼，測試一下

首先引入頭文件，然后定義一些用到的量

/* Private includes ----------------------------------------------------------*/
/* USER CODE BEGIN Includes */
#include "arm_math.h"
#include "arm_const_structs.h"
/* USER CODE END Includes */

/* Private typedef -----------------------------------------------------------*/
/* USER CODE BEGIN PTD */
#define Half_FFT_LENGTH	64
#define FFT_LENGTH      128				// FFT長度，默認是1024點FFT


uint16_t ADC_Value[FFT_LENGTH]; 		// ADC轉換結果
float fft_inputbuf[FFT_LENGTH*2];		// FFT輸入數組
float fft_outputbuf[FFT_LENGTH];		// FFT輸出數組
/* USER CODE END PTD */

然后編寫FFT變換的函數

void FFT(void)
{
	int i = 0;
	
	for(i=0;i < FFT_LENGTH;i++)
	{
		fft_inputbuf[i*2] = ADC_Value[i];
		fft_inputbuf[2*i +1] = 0;
	}
	
	arm_cfft_f32(&arm_cfft_sR_f32_len128,fft_inputbuf,0,1); // 執行FFT變換，arm_cfft_sR_f32_len128為宏，定義旋轉因子
	arm_cmplx_mag_f32(fft_inputbuf,fft_outputbuf,FFT_LENGTH);    // 把運算結果復數求模得幅值
}

這里要說明一下

DSP庫提供的fft變換，是可以同時進行正反變換的，為了兼容反變換，我們在進行FFT時，要將虛部賦值為0，即fft_inputbuf的內容應該為：ADC_Value1, 0, ADC_Value2, 0 .... 有效的值只有fft_inputbuf長度的一半，要進行的FFT也只有fft_inputbuf的一半，所以，要進行1024點的fft變換，就需要輸入2048長度的數組。

如果進行的是1024點的FFT變換，則得到1024長度的結果數組。

總諧波失真（THD）計算

首先，我們要知道經過FFT變換后得到的數組有什么物理意義，每個數組的值代表什么，這里直接說結論

\[設f_s 為采樣頻率，N為FFT點數（進行變換的FFT序列的長度），進行FFT變換后得數組下標為k的值代表的頻率為f_k 則： \]

\[f_k = \frac{k}{N}*f_s \]

也就是說，如果ADC的采樣頻率為1024Hz 進行2048點的FFT變換，則FFT變換后數組下標為5的位置對應的頻率為$$ f_5 = \frac{5}{1024} *2048 = 10hz$$ 也就是說，fft_outputbuf[5] 值的大小代表信號中10hz信號強度的大小。

明白了這個道理后，再來看如何計算總諧波失真，根據公式：

\[ THD = \frac{\sqrt{\sum_{n=0}^{\infty}{G_{n}^{2}}}}{G_0} (G_0為基波，G_n 為高次諧波) \]

我們需要知道基波和諧波分量的大小，就可以輕易的計算出來THD的大小

對E題來說，信號的頻率是固定的，只要采樣頻率\(f_s\)固定，點數N固定就可以很方便的得到基波分量，和諧波分量的位置

可以通過定時器來產生固定頻率的采樣信號，由於原信號為1k，諧波分量為2k、3k、4k、5k、6k... 采樣頻率應該大於最高頻率的3-5倍，

但，本題中我在題目中加入了波形顯示，根據采樣算出來的30k的采樣信號在LCD上的顯示效果不好，同時又為了方便計算，我的采樣頻率選擇為102.4kHz, 通過計算可以知道，fft_outputbuf[10]的大小就是基波分量的大小，二次諧波在fft_outputbuf[20]、三次諧波在fft_ouptubuf[30]... 依次類推， 這樣總諧波失真就可以很簡單的計算出來了。

void THD(void)
{
	thd_basic = fft_outputbuf[10];

	u[0] = fft_outputbuf[20];
	u[1] = fft_outputbuf[30];
	u[2] = fft_outputbuf[40];
	u[3] = fft_outputbuf[50];
	u[4] = fft_outputbuf[60];

	arm_power_f32(&u[0], 4, &sum);
	arm_sqrt_f32(sum, &thd_high);

	THD = thd_high / thd_basic;

}