完整版教程下載地址：http://www.armbbs.cn/forum.php?mod=viewthread&tid=94547
第29章 STM32F407移植匯編定點FFT庫（64點，256點和1024點）
本章主要講解ST官方匯編FFT庫的應用，包括1024點，256點和64點FFT的實作，
目錄
29.1 匯編FFT庫說明
29.1.1 描述
29.1.2 匯編庫的移植
29.2 函式cr4_fft_1024_stm32的使用(含幅頻和相頻回應)
29.3 函式cr4_fft_256_stm32的使用
29.4 函式cr4_fft_64_stm32的使用
29.5 實驗例程說明（MDK）
29.6 實驗例程說明（IAR）
29.7 總結
29.1 匯編FFT庫說明
29.1.1 描述
這個匯編的FFT庫是來自STM32F10x DSP library，由于是匯編實作的，而且是基4演算法，所以實作FFT在速度上比較快，
如果x[N]是采樣信號的話，使用FFT時必須滿足如下兩條：
- N得滿足4n（n =1，2, 3…..），也就是以4為基數，
- 采樣信號必須是32位資料，高16位存實部，低16位存虛部（這個是針對大端模式），小端模式是高位存虛部，低位存虛部，一般常用的是小端模式，
匯編FFT的實作主要包括以下三個函式：
- cr4_fft_64_stm32 ：實作64點FFT，
- cr4_fft_256_stm32 ：實作256點FFT，
- cr4_fft_1024_stm32 ： 實作1024點FFT，
29.1.2 匯編庫的移植
注：這里以MDK為例進行說明，IAR是一樣的，
這個匯編庫的移植比較簡單，從本章配套例子User檔案夾復制fft檔案夾到自己的工程：
注意路徑\User\fft\src\asm下有三個檔案夾，分布是arm，gcc和iar，其中arm可用于MDK，gcc可用于Embedded Studio，iar可用于IAR FOR ARM，
三個檔案夾里面都是如下幾個檔案，只是用于不用的編譯器：
然后把FFT源檔案的三個FFT匯編檔案和兩個頭檔案添加上即可，添加后效果如下(注意不同編譯器添相應匯編檔案)：
相應檔案添加后還有最重要一條，要把stm32_dsp.h檔案中的STM32F4的頭檔案：
最后別忘了添加路徑：
經過上面的操作，匯編FFT庫的移植就完成了，
29.2 函式cr4_fft_1024_stm32的使用(含幅頻和相頻回應)
cr4_fft_1024_stm32用于實作1024點資料的FFT計算，下面通過在開發板上運行這個函式并計算幅頻相應，然后再與Matlab計算的結果做對比，
uint32_t input[1024], output[1024], Mag[1024];/* 輸入，輸出和幅值 */
float32_t Phase[1024]; /* 相位*/
/*
*********************************************************************************************************
* 函 數 名: PowerMag
* 功能說明: 求模值
* 形 參：_usFFTPoints FFT點數
* 返 回 值: 無
*********************************************************************************************************
*/
void PowerMag(uint16_t _usFFTPoints)
{
int16_t lX,lY;
uint16_t i;
float32_t mag;
/* 計算幅值 */
for (i=0; i < _usFFTPoints; i++)
{
lX= (output[i]<<16)>>16; /* 實部*/
lY= (output[i]>> 16); /* 虛部 */
arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
Mag[i]= mag*2; /* 求模后乘以2才是實際模值，直流分量不需要乘2 */
}
/* 由于上面多乘了2，所以這里直流分量要除以2 */
Mag[0] = Mag[0]>>1;
}
/*
*********************************************************************************************************
* 函 數 名: Power_Phase_Radians
* 功能說明: 求相位
* 形 參：_usFFTPoints FFT點數， uiCmpValue 閥值
* 返 回 值: 無
*********************************************************************************************************
*/
void Power_Phase_Radians(uint16_t _usFFTPoints, uint32_t _uiCmpValue)
{
int16_t lX, lY;
uint16_t i;
float32_t phase;
float32_t mag;
for (i=0; i <_usFFTPoints; i++)
{
lX= (output[i]<<16)>>16; /* 實部 */
lY= (output[i] >> 16); /* 虛部 */
phase = atan2(lY, lX); /* atan2求解的結果范圍是(-pi, pi], 弧度制 */
arm_sqrt_f32((float32_t)(lX*lX+ lY*lY), &mag); /* 求模 */
if(_uiCmpValue > mag)
{
Phase[i] = 0;
}
else
{
Phase[i] = phase* 180.0f/3.1415926f; /* 將求解的結果由弧度轉換為角度 */
}
}
}
/*
*********************************************************************************************************
* 函 數 名: DSP_FFTPhase
* 功能說明: 1024點FFT的相位求解
* 形 參：無
* 返 回 值: 無
*********************************************************************************************************
*/
void DSP_FFTPhase(void)
{
uint16_t i;
/* 獲得1024個采樣點 */
for (i = 0; i < 1024; i++)
{
/* 波形是由直流分量，50Hz正弦波組成，波形采樣率1024 */
input[i] = 1024 + 1024*cos(2*3.1415926f*50*i/1024 + 3.1415926f/3);
}
/* 計算1024點FFT
output：輸出結果，高16位是虛部，低16位是實部，
input ：輸入資料，高16位是虛部，低16位是實部，
第三個引數必須是1024，
*/
cr4_fft_1024_stm32(output, input, 1024);
/* 求幅值 */
PowerMag(1024);
/* 列印輸出結果 */
for (i = 0; i < 1024; i++)
{
printf("%d\r\n", Mag[i]);
}
printf("=========================================\r\n");
/* 求相頻 */
Power_Phase_Radians(1024, 100);
/* 列印輸出結果 */
for (i = 0; i < 1024; i++)
{
printf("%f\r\n", Phase[i]);
}
}運行函式DSP_FFTPhase可以通過串口列印出計算的模值和相角，下面我們就通過Matlab計算的模值和相角跟cr4_fft_1024_stm32計算的做對比，
對比前需要先將串口列印出的資料加載到Matlab中，并給這個陣列起名sampledata，加載方法在前面的教程的第13章13.6小結已經講解，這里不做贅述了，Matlab中運行的代碼如下：
Fs = 1024; % 采樣率
N = 1024; % 采樣點數
n = 0:N-1; % 采樣序列
t = 0:1/Fs:1-1/Fs; % 時間序列
f = n * Fs / N; %真實的頻率
%波形是由直流分量，50Hz正弦波正弦波組成
x = 1024 + 1024*cos(2*pi*50*t + pi/3) ;
y = fft(x, N); %對原始信號做FFT變換
Mag = abs(y);
subplot(2,2,1);
MagAct = Mag *2 / N;
MagAct(1) = MagAct(1)/2;
plot(f, MagAct);
title('Matlab計算幅頻回應');
xlabel('頻率');
ylabel('賦值');
subplot(2,2,2);
realvalue = real(y);
imagvalue = imag(y);
plot(f, atan2(imagvalue, realvalue)*180/pi.*(Mag>=1024*20));
title('Matlab計算相頻回應');
xlabel('頻率');
ylabel('相角');
subplot(2,2,3);
plot(f, sampledata1); %繪制STM32計算的幅頻相應
title('STM32計算幅頻回應');
xlabel('頻率');
ylabel('賦值');
subplot(2,2,4);
plot(f, sampledata2); %繪制STM32計算的相頻相應
title('STM32計算相頻回應');
xlabel('頻率');
ylabel('相角');運行Matlab后的輸出結果如下：
從上面的對比結果中可以看出，Matlab和函式cr4_fft_1024_stm32計算的結果基本是一直的，幅頻回應求出的幅值和相頻回應中的求出的初始相角都是沒問題的，
29.3 函式cr4_fft_256_stm32的使用
cr4_fft_256_stm32和cr4_fft_1024_stm32的用法是一樣的，下面通過一個實體進行說明：
/*
*********************************************************************************************************
* 函 數 名: DSP_FFT256
* 功能說明: 256點FFT實作
* 形 參：無
* 返 回 值: 無
*********************************************************************************************************
*/
void DSP_FFT256(void)
{
uint16_t i;
/* 獲得256個采樣點 */
for (i = 0; i < 256; i++)
{
/* 波形是由直流分量，50Hz正弦波和20Hz正弦波組成，波形采樣率200Hz */
input[i] = 1024 + 1024*sin(2*3.1415926f*50*i/200) + 512*sin(2*3.1415926f*20*i/200) ;
}
/* 計算256點FFT
output：輸出結果，高16位是虛部，低16位是實部，
input ：輸入資料，高16位是虛部，低16位是實部，
第三個引數必須是1024，
*/
cr4_fft_256_stm32(output, input, 256);
/* 求幅值 */
PowerMag(256);
/* 列印輸出結果 */
for (i = 0; i < 256; i++)
{
printf("%d\r\n", Mag[i]);
}
}運行函式DSP_FFT256可以通過串口列印出計算的模值，下面我們就通過Matlab計算的模值跟cr4_fft_256_stm32計算的模值做對比，
對比前需要先將串口列印出的資料加載到Matlab中，并給這個陣列起名sampledata， Matlab中運行的代碼如下：
Fs = 200; % 采樣率
N = 256; % 采樣點數
n = 0:N-1; % 采樣序列
t = 0:1/Fs:1-1/Fs; % 時間序列
f = n * Fs / N; %真實的頻率
%波形是由直流分量，50Hz正弦波和20Hz正弦波組成
x = 1024 + 1024*sin(2*pi*50*t) + 512*sin(2*pi*20*t) ;
y = fft(x, N); %對原始信號做FFT變換
subplot(2,1,1);
Mag = abs(y);
MagAct = Mag *2 / N;
MagAct(1) = MagAct(1)/2;
plot(f, MagAct); %繪制幅頻相應曲線
title('Matlab計算結果');
xlabel('頻率');
ylabel('幅度');
subplot(2,1,2);
plot(f, sampledata); %繪制STM32計算的幅頻相應
title('STM32計算結果');
xlabel('頻率');
ylabel('幅度');運行Matlab后的輸出結果如下：
從上面的對比結果中可以看出，Matlab和函式cr4_fft_256_stm32計算的結果基本是一直的，但頻率泄露略多，
29.4 函式cr4_fft_64_stm32的使用
cr4_fft_64_stm32和cr4_fft_1024_stm32的用法也是一樣的，下面通過一個實體進行說明：
/*
*********************************************************************************************************
* 函 數 名: DSP_FFT64
* 功能說明: 64點FFT實作
* 形 參：無
* 返 回 值: 無
*********************************************************************************************************
*/
void DSP_FFT64(void)
{
uint16_t i;
/* 獲得64個采樣點 */
for (i = 0; i < 64; i++)
{
/* 波形是由直流分量，5Hz正弦波和10Hz正弦波組成，波形采樣率60Hz */
input[i] = 1024 + 1024*sin(2*3.1415926f*5*i/60) + 512*sin(2*3.1415926f*10*i/60) ;
}
/* 計算64點FFT
output：輸出結果，高16位是虛部，低16位是實部，
input ：輸入資料，高16位是虛部，低16位是實部，
第三個引數必須是1024，
*/
cr4_fft_64_stm32(output, input, 64);
/* 求幅值 */
PowerMag(64);
/* 列印輸出結果 */
for (i = 0; i < 64; i++)
{
printf("%d\r\n", Mag[i]);
}
}運行函式DSP_FFT64可以通過串口列印出計算的模值，下面我們就通過Matlab計算的模值跟cr4_fft_64_stm32計算的模值做對比，
對比前需要先將串口列印出的資料加載到Matlab中，并給這個陣列起名sampledata，Matlab中運行的代碼如下：
Fs = 60; % 采樣率
N = 64; % 采樣點數
n = 0:N-1; % 采樣序列
t = 0:1/Fs:1-1/Fs; % 時間序列
f = n * Fs / N; % 真實的頻率
%波形是由直流分量，5Hz正弦波和10Hz正弦波組成
x = 1024 + 1024*sin(2*pi*5*t) + 512*sin(2*pi*10*t) ;
y = fft(x, N); %對原始信號做FFT變換
subplot(2,1,1);
Mag = abs(y);
MagAct = Mag *2 / N;
MagAct(1) = MagAct(1)/2;
plot(f, MagAct); %繪制幅頻相應曲線
title('Matlab計算結果');
xlabel('頻率');
ylabel('幅度');
subplot(2,1,2);
plot(f, sampledata); %繪制STM32計算的幅頻相應
title('STM32計算結果');
xlabel('頻率');
ylabel('幅度');運行Matlab后的輸出結果如下：
從上面的對比結果中可以看出，Matlab和函式cr4_fft_64_stm32計算的結果基本是一直的，但是計算的效果都比較差，主要是因為采樣點數太少，
29.5 實驗例程說明（MDK）
配套例子：
V5-219_STM32F407移植ST匯編定點FFT庫(64點，256點和1024點)
實驗目的：
- 學習ST匯編定點FFT庫（64點，256點和1024點）
實驗內容：
- 啟動一個自動重裝軟體定時器，每100ms翻轉一次LED2，
- 按下按鍵K1，串口列印1024點FFT的幅頻回應和相頻回應，
- 按下按鍵K2，串口列印256點FFT的幅頻回應，
- 按下按鍵K3，串口列印64點FFT的幅頻回應，
使用AC6注意事項
特別注意附件章節C的問題
上電后串口列印的資訊：
波特率 115200，資料位 8，奇偶校驗位無，停止位 1，
RTT方式列印資訊：
程式設計：
系統堆疊大小分配：
硬體外設初始化
硬體外設的初始化是在 bsp.c 檔案實作：
/*
*********************************************************************************************************
* 函 數 名: bsp_Init
* 功能說明: 初始化所有的硬體設備，該函式配置CPU暫存器和外設的暫存器并初始化一些全域變數，只需要呼叫一次
* 形 參：無
* 返 回 值: 無
*********************************************************************************************************
*/
void bsp_Init(void)
{
/*
STM32F407 HAL 庫初始化，此時系統用的還是F407自帶的16MHz，HSI時鐘:
- 呼叫函式HAL_InitTick，初始化滴答時鐘中斷1ms，
- 設定NVIV優先級分組為4，
*/
HAL_Init();
/*
配置系統時鐘到168MHz
- 切換使用HSE，
- 此函式會更新全域變數SystemCoreClock，并重新配置HAL_InitTick，
*/
SystemClock_Config();
/*
Event Recorder：
- 可用于代碼執行時間測量，MDK5.25及其以上版本才支持，IAR不支持，
- 默認不開啟，如果要使能此選項，務必看V5開發板用戶手冊第8章
*/
#if Enable_EventRecorder == 1
/* 初始化EventRecorder并開啟 */
EventRecorderInitialize(EventRecordAll, 1U);
EventRecorderStart();
#endif
bsp_InitKey(); /* 按鍵初始化，要放在滴答定時器之前，因為按鈕檢測是通過滴答定時器掃描 */
bsp_InitTimer(); /* 初始化滴答定時器 */
bsp_InitUart(); /* 初始化串口 */
bsp_InitExtIO(); /* 初始化擴展IO */
bsp_InitLed(); /* 初始化LED */
}主功能：
主程式實作如下操作：
- 啟動一個自動重裝軟體定時器，每100ms翻轉一次LED2，
- 按下按鍵K1，串口列印1024點FFT的幅頻回應和相頻回應，
- 按下按鍵K2，串口列印256點FFT的幅頻回應，
- 按下按鍵K3，串口列印64點FFT的幅頻回應，
/*
*********************************************************************************************************
* 函 數 名: main
* 功能說明: c程式入口
* 形 參: 無
* 返 回 值: 錯誤代碼(無需處理)
*********************************************************************************************************
*/
int main(void)
{
uint8_t ucKeyCode; /* 按鍵代碼 */
bsp_Init(); /* 硬體初始化 */
PrintfLogo(); /* 列印例程資訊到串口1 */
PrintfHelp(); /* 列印操作提示資訊 */
bsp_StartAutoTimer(0, 100); /* 啟動1個100ms的自動重裝的定時器 */
/* 進入主程式回圈體 */
while (1)
{
bsp_Idle(); /* 這個函式在bsp.c檔案，用戶可以修改這個函式實作CPU休眠和喂狗 */
if (bsp_CheckTimer(0)) /* 判斷定時器超時時間 */
{
/* 每隔500ms 進來一次 */
bsp_LedToggle(4); /* 翻轉LED2的狀態 */
}
ucKeyCode = bsp_GetKey(); /* 讀取鍵值, 無鍵按下時回傳 KEY_NONE = 0 */
if (ucKeyCode != KEY_NONE)
{
switch (ucKeyCode)
{
case KEY_DOWN_K1: /* K1鍵按下 */
DSP_FFTPhase();
break;
case KEY_DOWN_K2: /* K2鍵按下 */
DSP_FFT256();
break;
case KEY_DOWN_K3: /* K3鍵按下 */
DSP_FFT64();
break;
default:
/* 其它的鍵值不處理 */
break;
}
}
}
}29.6 實驗例程說明（IAR）
配套例子：
V5-219_STM32F407移植ST匯編定點FFT庫(64點，256點和1024點)
實驗目的：
- 學習ST匯編定點FFT庫（64點，256點和1024點）
實驗內容：
- 啟動一個自動重裝軟體定時器，每100ms翻轉一次LED2，
- 按下按鍵K1，串口列印1024點FFT的幅頻回應和相頻回應，
- 按下按鍵K2，串口列印256點FFT的幅頻回應，
- 按下按鍵K3，串口列印64點FFT的幅頻回應，
上電后串口列印的資訊：
波特率 115200，資料位 8，奇偶校驗位無，停止位 1，
RTT方式列印資訊：
程式設計：
系統堆疊大小分配：
硬體外設初始化
硬體外設的初始化是在 bsp.c 檔案實作：
/*
*********************************************************************************************************
* 函 數 名: bsp_Init
* 功能說明: 初始化所有的硬體設備，該函式配置CPU暫存器和外設的暫存器并初始化一些全域變數，只需要呼叫一次
* 形 參：無
* 返 回 值: 無
*********************************************************************************************************
*/
void bsp_Init(void)
{
/*
STM32F407 HAL 庫初始化，此時系統用的還是F407自帶的16MHz，HSI時鐘:
- 呼叫函式HAL_InitTick，初始化滴答時鐘中斷1ms，
- 設定NVIV優先級分組為4，
*/
HAL_Init();
/*
配置系統時鐘到168MHz
- 切換使用HSE，
- 此函式會更新全域變數SystemCoreClock，并重新配置HAL_InitTick，
*/
SystemClock_Config();
/*
Event Recorder：
- 可用于代碼執行時間測量，MDK5.25及其以上版本才支持，IAR不支持，
- 默認不開啟，如果要使能此選項，務必看V5開發板用戶手冊第8章
*/
#if Enable_EventRecorder == 1
/* 初始化EventRecorder并開啟 */
EventRecorderInitialize(EventRecordAll, 1U);
EventRecorderStart();
#endif
bsp_InitKey(); /* 按鍵初始化，要放在滴答定時器之前，因為按鈕檢測是通過滴答定時器掃描 */
bsp_InitTimer(); /* 初始化滴答定時器 */
bsp_InitUart(); /* 初始化串口 */
bsp_InitExtIO(); /* 初始化擴展IO */
bsp_InitLed(); /* 初始化LED */
}主功能：
主程式實作如下操作：
- 啟動一個自動重裝軟體定時器，每100ms翻轉一次LED2，
- 按下按鍵K1，串口列印1024點FFT的幅頻回應和相頻回應，
- 按下按鍵K2，串口列印256點FFT的幅頻回應，
- 按下按鍵K3，串口列印64點FFT的幅頻回應，
/*
*********************************************************************************************************
* 函 數 名: main
* 功能說明: c程式入口
* 形 參: 無
* 返 回 值: 錯誤代碼(無需處理)
*********************************************************************************************************
*/
int main(void)
{
uint8_t ucKeyCode; /* 按鍵代碼 */
bsp_Init(); /* 硬體初始化 */
PrintfLogo(); /* 列印例程資訊到串口1 */
PrintfHelp(); /* 列印操作提示資訊 */
bsp_StartAutoTimer(0, 100); /* 啟動1個100ms的自動重裝的定時器 */
/* 進入主程式回圈體 */
while (1)
{
bsp_Idle(); /* 這個函式在bsp.c檔案，用戶可以修改這個函式實作CPU休眠和喂狗 */
if (bsp_CheckTimer(0)) /* 判斷定時器超時時間 */
{
/* 每隔500ms 進來一次 */
bsp_LedToggle(4); /* 翻轉LED2的狀態 */
}
ucKeyCode = bsp_GetKey(); /* 讀取鍵值, 無鍵按下時回傳 KEY_NONE = 0 */
if (ucKeyCode != KEY_NONE)
{
switch (ucKeyCode)
{
case KEY_DOWN_K1: /* K1鍵按下 */
DSP_FFTPhase();
break;
case KEY_DOWN_K2: /* K2鍵按下 */
DSP_FFT256();
break;
case KEY_DOWN_K3: /* K3鍵按下 */
DSP_FFT64();
break;
default:
/* 其它的鍵值不處理 */
break;
}
}
}
}29.7 總結
本章節主要講解了匯編FFT的1024點，256點和64點使用方法，有興趣的可以深入了解匯編代碼的實作，
轉載請註明出處，本文鏈接：https://www.uj5u.com/qita/286999.html
標籤：其他
上一篇：stm32f103 FatFs檔案系統移植野火的代碼到正點原子開發板上
下一篇：RS-485通信協議簡介