關于adc采樣與fft轉換之間的關系與應用

寫在前面

最近在要實作一個功能，對外部波形進行采樣，然後用官方fft變換對外部波形進行諧波分析。但是本人數學不好，不喜歡研究算法，不喜歡看一堆公式，不想要知道原理，隻想知道如何應用。在網上找了很久都沒有找到想要的demo，大部分都解釋得不夠全面，自己研究後，上傳一包源碼給大家分享。本人能力有限，關于後文對其中的解釋，有不正确的地方，還請大家指正。

如果有對傅裡葉變換原理感興趣的朋友，強烈推薦大家看這一位大佬的部落格，對了解傅裡葉分析的意義很有幫助：https://www.cnblogs.com/h2zZhou/p/8405717.html

附上完整的工程連結，需要的朋友請自行下載下傳（代碼中的注釋部分，采用ANSI編碼格式）：https://github.com/WangYooNestik/STM32F10X-ADC-FFT

直接進入主題

1.說一下總的架構

1.目的：用ADC采樣值進行fft運算，求取采樣後波形相關屬性（本例主要分析50Hz正弦波，及高/低次諧波）
2.方案：計劃ADC采樣頻率為512Hz，采樣1024個點（也就是采樣2秒鐘）進行一次fft運算，使頻域分辨率達到0.5Hz
3.配置：由于我用的系統時鐘是72MHz，是以ADC本身不能配置到512Hz，采用tim3定時觸發采樣（當系統時鐘為72MHz時，ADC時鐘最慢是9MHz，采樣周期最大約為256個時鐘周期，是以ADC本身采樣頻率最慢約為35156.25Hz（大于我需要的512Hz），是以選擇TIM定時采樣）
4.其他：
- 為什麼要選512Hz：為了滿足采樣定理（fs > 100Hz）；又因為我想要分析高次諧波；又因為我選了1024的fft變換，為了使頻域分辨率識别友善（0.5Hz）。這句話有點繞，後文有詳細解釋。
- 本來計劃采樣5120個點（使頻域分辨率達到0.1Hz），但是單片機RAM空間不夠，有條件的可以嘗試
5.相關知識（很重要，本文用到的知識點都在這）：
- 傅裡葉變換的目的：求取被測波形的幅頻特性/相頻特性
- 采樣定理：采樣頻率需大于被采樣波形頻率的兩倍
- 頻域分辨率 = fs/N，（fs：采樣頻率，N：采樣點數）
- 幅頻特性分析範圍（機關：Hz）：0 ~ fs/2，（fs：采樣頻率）
- 幅頻特性序列：FFT_Mag[i] 代表的意思，頻率分量為 (i*(fs/N)) 的幅值為 FFT_Mag[i]

2.說一下采樣頻率及點數的選擇（重點）

由于fft變換需要考慮多個變量對最終結果的影響，比較難辦；我建議，對采樣頻率及點數的選擇按以下步驟進行分析：

1.采樣點數分析：
- 不考慮耗時，采樣點數越多越好。
- 看你的平台的資源情況，分析整個項目的RAM大小及占用情況，先選擇你準備用多少點進行fft變。當你的fs确定下來，采樣點數越多，你的頻域分辨率越小，你最終的幅頻特性序列就越精确。當然我們一般認為越精确越好。（舉個例子：我分析了1024個點，則fft變換輸入數組大小 = 1024 * 4Byte = 4K，輸出數組大小 = 1024 * 4Byte = 4K，這樣就已經占用8KRAM）
2.采樣頻率分析：
- 滿足采樣定理：采樣頻率需大于被采樣波形頻率的兩倍（假設你需要采樣的波形頻率為50Hz，則采樣頻率必須在100Hz以上）。
- 如果你隻想分析一個特定頻率的波形（比如50Hz），不考慮它的高（>50Hz）/低（<50Hz）次諧波，隻需要滿足采樣定理即可；
- 如果你需要分析特定頻率及低次諧波，隻需要滿足采樣定理即可；
- 如果你需要分析特定頻率，及它的高/低次諧波。那你需要優先考慮你的采樣頻率對你需要分析的最高次諧波，要滿足采樣定理；也就是fs要 > 2 * 最高次諧波頻率。
  
  反過來說，就是要考慮你的幅頻特性分析範圍。
3.采樣精度确認：
- 确定了以上兩點後，計算頻率分辨率（頻域分辨率 = fs/N），看是否滿足你的精度要求。

3.官方給出了fft庫函數啥意思，咋用？（官方庫在哪裡下載下傳？我也不知道）

/* 64 points*/
void cr4_fft_64_stm32(void *pssOUT, void *pssIN, uint16_t Nbin);
/* 256 points */
void cr4_fft_256_stm32(void *pssOUT, void *pssIN, uint16_t Nbin);
/* 1024 points */
void cr4_fft_1024_stm32(void *pssOUT, void *pssIN, uint16_t Nbin);

官方給出了3個fft函數，分别是64、256、1024個點的fft變換（就是采樣的樣本是64、256、1024個樣本）。

pssIN：是adc采樣的樣本

pssOUT：調用官方函數後，輸出的傅裡葉序列

Nbin：樣本的數量（采了多少個點）

為什麼是64、256、1024這個3個值？因為官方使用的是基4蝶形算法（啥意思？），簡單了解要使用官方fft庫函數，必須遵循采樣點數Nbin，是4的n次方。

4.傅裡葉變換後，輸出了個啥，咋用？

我們需要有個基礎知識，傅裡葉變換的目的：求取幅頻特性/相頻特性

fft變換後，輸出的是一個傅裡葉序列（怎麼用？）。傅裡葉序列本身，不是我們能夠肉眼分析的東西，我們還需要對傅裡葉序列進行計算，求取幅頻特性/相頻特性序列。（下面這段代碼就是求取幅頻特性）

#define SAMPLS_NUM 1024
u32 FFT_OutData[SAMPLS_NUM] = {0};		//fft輸出序列
u32 FFT_Mag[SAMPLS_NUM/2] = {0};		//幅頻特性序列（序号代表頻率分量，值代表幅值大小。由于FFT的頻譜結果是關于奈奎斯特頻率對稱的，是以隻計算一半的點即可）
/********************************************************************************************************
函數名稱:GetPowerMag()
函數功能:計算各次諧波幅值
參數說明:
備　　注:先将ADC_FFT_OutData分解成實部(X)和虛部(Y)，然後計算幅頻特性序列FFT_Mag
         本函數參考網頁：https://wenku.baidu.com/view/08ccee0984868762cbaed532.html，關于幅頻特性計算部分
**********************************************************************************************************/
void GetPowerMag(void)
{
	signed short lX,lY;
	float X,Y,Mag;
	unsigned short i;

	for(i=0; i<SAMPLS_NUM/2; i++)
	{
		lX = (FFT_OutData[i] << 16) >> 16;
		lY = (FFT_OutData[i] >> 16);
		
		X = SAMPLS_NUM * ((float)lX) / 32768;
		Y = SAMPLS_NUM * ((float)lY) / 32768;
		
		Mag = sqrt(X * X + Y * Y) / SAMPLS_NUM;
		
		if(i == 0)
			FFT_Mag[i] = (unsigned long)(Mag * 32768);
		else
			FFT_Mag[i] = (unsigned long)(Mag * 65536);
		//printf("%ld\r\n",lBufMagArray[i]);
	}
	//printf("\r\n\r\n");
}

FFT_Mag[SAMPLS_NUM/2]，就是幅頻特性序列，到此傅裡葉變換的目的就達到了。

我輸入的是50Hz正弦波，可以看到FFT_Mag[100]這個值，遠大于其他值。證明頻率分量為50Hz（頻域分辨率（0.5Hz） * 100）的諧波，幅值最大，是以可以知道fft後是正确的。

5.最後貼上源碼

這個玩意兒好像下載下傳源碼要積分，我就直接貼在後面了，講得不是很清楚大家見諒。

ADC配置相關

#define ADC_CHANNEL_NUMS 	3
#define SAMPLS_NUM        	1024

#define TIM2_PERIOD 1953
#define TIM3_PERIOD 1953


u16 ADC_SourceData[SAMPLS_NUM][ADC_CHANNEL_NUMS] = {0};



void ADC_GPIO_Configuration(void)
{
	GPIO_InitTypeDef GPIO_InitStructure;

	RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA, ENABLE);				//使能GPIOA時鐘

	GPIO_InitStructure.GPIO_Pin = GPIO_Pin_1|GPIO_Pin_4|GPIO_Pin_5;		//管腳1/4/5
	GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AIN;						//模拟輸入模式
	GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;
	GPIO_Init(GPIOA, &GPIO_InitStructure);								//GPIO組
}

void ADC_TIM2_GPIO_Configuration(void)        //ADC配置函數
{
	GPIO_InitTypeDef GPIO_InitStructure;

	RCC_APB2PeriphClockCmd(RCC_APB2Periph_GPIOA, ENABLE);   //使能GPIOA時鐘

	GPIO_InitStructure.GPIO_Pin = GPIO_Pin_2;				//TIM2_CH3
	GPIO_InitStructure.GPIO_Mode = GPIO_Mode_AF_PP;
	GPIO_InitStructure.GPIO_Speed = GPIO_Speed_50MHz;
	GPIO_Init(GPIOA, &GPIO_InitStructure);
}

void ADC_TIM2_Configuration(void)
{ 
	TIM_TimeBaseInitTypeDef TIM_TimeBaseStructure;
	TIM_OCInitTypeDef TIM_OCInitStructure;

	ADC_TIM2_GPIO_Configuration();
	
	RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM2,ENABLE);

	TIM_DeInit(TIM2);

	//配置tim2觸發頻率為512Hz
	TIM_TimeBaseStructure.TIM_Period = TIM2_PERIOD - 1;				//設定2ms一次TIM2比較的周期
	TIM_TimeBaseStructure.TIM_Prescaler = 71;						//系統主頻72M，這裡分頻71，相當于1000K的定時器2時鐘
	TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV1;
	TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;
	TIM_TimeBaseInit(TIM2, & TIM_TimeBaseStructure);

	TIM_OCInitStructure.TIM_OCMode = TIM_OCMode_PWM1;
	TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;	//TIM_OutputState_Disable;
	TIM_OCInitStructure.TIM_Pulse = (TIM2_PERIOD - 1) / 2;
	TIM_OCInitStructure.TIM_OCPolarity = TIM_OCPolarity_Low;		//如果是PWM1要為Low，PWM2則為High
	TIM_OC3Init(TIM2, & TIM_OCInitStructure);
	
	//TIM_InternalClockConfig(TIM2);
	
	//TIM_OC3PreloadConfig(TIM2,TIM_OCPreload_Enable);
	
	//TIM_UpdateDisableConfig(TIM2,DISABLE);

	//TIM_Cmd(TIM2, ENABLE);//最後面打開定時器使能

	//TIM_CtrlPWMOutputs(TIM2,ENABLE);
}

void ADC_TIM3_Configuration(void)
{
	TIM_TimeBaseInitTypeDef TIM_TimeBaseInitStructure; 
	
	RCC_APB1PeriphClockCmd(RCC_APB1Periph_TIM3,ENABLE);

	TIM_TimeBaseInitStructure.TIM_Period = TIM3_PERIOD - 1;		//設定ADC_CHANNEL_NUMS * 512Hz采樣頻率
	TIM_TimeBaseInitStructure.TIM_Prescaler = 71;
	TIM_TimeBaseInitStructure.TIM_ClockDivision = TIM_CKD_DIV1;			//時鐘分頻，TIM3是通用定時器，基本定時器不用設定
	TIM_TimeBaseInitStructure.TIM_CounterMode = TIM_CounterMode_Up;		//向上掃描
	TIM_TimeBaseInit(TIM3, &TIM_TimeBaseInitStructure);
	
	TIM_SelectOutputTrigger(TIM3, TIM_TRGOSource_Update);				//選擇TRGO作為觸發源為定時器更新時間
	
	TIM_Cmd(TIM3,ENABLE);												//開啟定時器3
}

void ADC_DMA_NVIC_Configuration(void)
{
	NVIC_InitTypeDef NVIC_InitStructure;

	NVIC_InitStructure.NVIC_IRQChannel					 = DMA1_Channel1_IRQn;
	NVIC_InitStructure.NVIC_IRQChannelPreemptionPriority = 2;
	NVIC_InitStructure.NVIC_IRQChannelSubPriority		 = 1;
	NVIC_InitStructure.NVIC_IRQChannelCmd                = ENABLE;
	NVIC_Init(&NVIC_InitStructure);

	DMA_ClearITPendingBit(DMA1_IT_TC1);

	DMA_ITConfig(DMA1_Channel1,DMA1_IT_TC1,ENABLE);
}

void ADC_DMA_Configuration(void)
{
	DMA_InitTypeDef  DMA_InitStructure;

	RCC_AHBPeriphClockCmd(RCC_AHBPeriph_DMA1,ENABLE);

	DMA_DeInit(DMA1_Channel1);
	DMA_InitStructure.DMA_PeripheralBaseAddr = (u32)&ADC1->DR;
	DMA_InitStructure.DMA_MemoryBaseAddr     = (u32)ADC_SourceData;
	DMA_InitStructure.DMA_DIR                = DMA_DIR_PeripheralSRC;
	DMA_InitStructure.DMA_BufferSize         = ADC_CHANNEL_NUMS*SAMPLS_NUM;
	DMA_InitStructure.DMA_PeripheralInc      = DMA_PeripheralInc_Disable;
	DMA_InitStructure.DMA_MemoryInc          = DMA_MemoryInc_Enable;
	DMA_InitStructure.DMA_PeripheralDataSize = DMA_PeripheralDataSize_HalfWord;
	DMA_InitStructure.DMA_MemoryDataSize     = DMA_MemoryDataSize_HalfWord;
	DMA_InitStructure.DMA_Mode               = DMA_Mode_Circular;
	DMA_InitStructure.DMA_Priority           = DMA_Priority_High;
	DMA_InitStructure.DMA_M2M                = DMA_M2M_Disable;
	DMA_Init(DMA1_Channel1, &DMA_InitStructure);

	DMA_Cmd(DMA1_Channel1, ENABLE);//使能DMA	

	ADC_DMA_NVIC_Configuration();
}

void ADC_Init_Configuration(void)//ADC配置函數
{
	ADC_InitTypeDef  ADC_InitStructure;

	RCC_APB2PeriphClockCmd(RCC_APB2Periph_ADC1, ENABLE);
	
	RCC_ADCCLKConfig(RCC_PCLK2_Div6);

	ADC_DeInit(ADC1);
	ADC_InitStructure.ADC_Mode               = ADC_Mode_Independent;
	ADC_InitStructure.ADC_ScanConvMode       = ENABLE;								//如果是單通道，關閉掃描模式
	ADC_InitStructure.ADC_ContinuousConvMode = DISABLE;
	ADC_InitStructure.ADC_ExternalTrigConv   = ADC_ExternalTrigConv_T3_TRGO;
	ADC_InitStructure.ADC_DataAlign          = ADC_DataAlign_Right;
	ADC_InitStructure.ADC_NbrOfChannel       = ADC_CHANNEL_NUMS;
	ADC_Init(ADC1, &ADC_InitStructure);

	//==========================================================================   
	ADC_RegularChannelConfig(ADC1, ADC_Channel_1,  1, ADC_SampleTime_239Cycles5);	//AI_VS_A1
	ADC_RegularChannelConfig(ADC1, ADC_Channel_4,  2, ADC_SampleTime_239Cycles5);	//AI_VS_B1
	ADC_RegularChannelConfig(ADC1, ADC_Channel_5,  3, ADC_SampleTime_239Cycles5);	//AI_VS_C1

	ADC_DMACmd(ADC1, ENABLE);
	
	ADC_Cmd(ADC1, ENABLE);
	
	ADC_ResetCalibration(ADC1);
	while(ADC_GetResetCalibrationStatus(ADC1));
	ADC_StartCalibration(ADC1);
	while(ADC_GetCalibrationStatus(ADC1));

	ADC_ExternalTrigConvCmd(ADC1, ENABLE);
	//ADC_SoftwareStartConvCmd(ADC1, ENABLE);
}



/***********************************************************************************************************************************************
*adc初始化配置
*1.目     的：用ADC采樣值進行fft運算，求取采樣後波形相關屬性（本例主要分析50Hz正選波，及高/低次諧波）
*2.方     案：計劃ADC采樣頻率為512Hz，采樣1024個點進行fft運算，使頻域分辨率達到0.5Hz
*3.配     置：由于ADC本身不能配置到512Hz，采用tim3定時觸發采樣
*4.其     他：當系統時鐘為72MHz時，ADC時鐘最慢是9MHz，采樣周期最大約為256個時鐘周期，是以ADC本身采樣頻率最慢約為35156.25Hz（大于我需要的512Hz）
*            為什麼要選512Hz：由于采用标準fft基4算法（采樣點數必須是4的倍數），為了使頻域分辨率為識别友善（0.5Hz），又為了滿足采樣定理
*            本來計劃采樣5120個點（使頻域分辨率達到0.1Hz），但是單片機RAM空間不夠，有條件的可以嘗試
*5.相關知識：采樣定理：采樣頻率需大于被采樣波形頻率的兩倍
*            頻域分辨率 = fs/N，（fs：采樣頻率，N：采樣點數）
*            幅頻特性序列：FFT_Mag[i]代表的意思，頻率分量為i*(fs/N)的幅值為FFT_Mag[i]
************************************************************************************************************************************************/
void Adc_Init(void)
{
	ADC_GPIO_Configuration();
	
	//ADC_TIM2_Configuration();

	ADC_TIM3_Configuration();

	ADC_DMA_Configuration();

	ADC_Init_Configuration();
}

//ADC_DMA中斷服務程式
void DMA1_Channel1_IRQHandler(void)
{
	if(DMA_GetITStatus(DMA1_IT_TC1) != RESET)
	{
		global.adc_finish_fg = true;
		DMA_ClearITPendingBit(DMA1_IT_TC1);
	}
}

fft函數

u32 FFT_SourceData[SAMPLS_NUM] = {0};	//fft輸入序列
u32 FFT_OutData[SAMPLS_NUM] = {0};		//fft輸出序列
u32 FFT_Mag[SAMPLS_NUM/2] = {0};		//幅頻特性序列（序号代表頻率分量，值代表幅值大小。由于FFT的頻譜結果是關于奈奎斯特頻率對稱的，是以隻計算一半的點即可）


#if 0
/****************************************************************************************************
函數名稱:InitBufInArray()
函數功能:模拟采樣資料，采樣資料中包含3種頻率正弦波(350Hz，8400Hz，18725Hz)
參數說明:
備     注:在lBufInArray數組中，每個資料的高16位存儲采樣資料的實部，低16位存儲采樣資料的虛部(總是為0)
****************************************************************************************************/

u16 lBufInArray[SAMPLS_NUM];

void InitBufInArray(void)
{
	unsigned short i;
	float fx;
	for(i=0; i<SAMPLS_NUM; i++)
	{
		fx = 1500 * sin(PI2 * i * 55.6 / Fs);
		lBufInArray[i] = ((signed short)fx)<<16;
	}
}
#endif

/********************************************************************************************************
函數名稱:GetPowerMag()
函數功能:計算各次諧波幅值
參數說明:
備　　注:先将ADC_FFT_OutData分解成實部(X)和虛部(Y)，然後計算幅頻特性序列FFT_Mag
         本函數參考網頁：https://wenku.baidu.com/view/08ccee0984868762cbaed532.html，關于幅頻特性計算部分
**********************************************************************************************************/
void GetPowerMag(void)
{
	signed short lX,lY;
	float X,Y,Mag;
	unsigned short i;

	for(i=0; i<SAMPLS_NUM/2; i++)
	{
		lX = (FFT_OutData[i] << 16) >> 16;
		lY = (FFT_OutData[i] >> 16);
		
		X = SAMPLS_NUM * ((float)lX) / 32768;
		Y = SAMPLS_NUM * ((float)lY) / 32768;
		
		Mag = sqrt(X * X + Y * Y) / SAMPLS_NUM;
		
		if(i == 0)
			FFT_Mag[i] = (unsigned long)(Mag * 32768);
		else
			FFT_Mag[i] = (unsigned long)(Mag * 65536);
		//printf("%ld\r\n",lBufMagArray[i]);
	}
	//printf("\r\n\r\n");
}

/************************************************
*由于在fft運算過程中不允許源資料更新
*将源資料拷貝到另一塊記憶體
*此步驟可以被替換為：進行fft運算時，關閉adc轉換
*************************************************/
void Get_FFT_Source_Data(EN_FFT_CHANNEL channel_idx)
{
	u16 i,j;

	for(i=0; i<SAMPLS_NUM; i++)
	{
		FFT_SourceData[i] = (u32)ADC_SourceData[i][channel_idx];
	}
}

void FFT_test(void)
{
	//InitBufInArray();
	Get_FFT_Source_Data(FFT_CHANNEL_1);
	cr4_fft_1024_stm32(FFT_OutData, FFT_SourceData, SAMPLS_NUM);
	GetPowerMag();
}

主函數

int main(void)
{
	Adc_Init();

	while(1)
	{
		if(global.adc_finish_fg)
		{
			global.adc_finish_fg = false;
			FFT_test();				//512Hz采樣頻率，采樣1024點，可以看到2s鐘周期執行FFT_test();
		}
	}
}

STM32F1X+ADC+FFT應用寫在前面

寫在前面

直接進入主題

繼續閱讀

Matlab中FFT函數的源代碼示例

LabVIEW FPGA PCIe開發講解-7.1節：FPGA PCIe/PXIe基礎知識和概念概述

13-STM32F1 DAC輸出

9-STM32F1-輸入捕獲

5-STM32F1 外部中斷

22-STM32F1-WIFI子產品ESP8266

6-STM32F1 -獨立看門狗

stm32 usart 單線半雙工序列槽 DMA發送最後一個位元組發不出來問題

DSP C2000系列TMS320F28335學習之ADC

POJ 2389 Bull Math （FFT）

STM32F1學習筆記(二)STM32學習之GPIO

[轉載]20世紀十大算法

STFT filter bankSTFT filter bank

初出茅廬的小李第16篇部落格之三路ADC采集資料DMA方式轉換

相位相關算法

STM32F103雙向可控矽調節220V調光簡介：原理：程式邏輯：注意：