可用於嵌入式計算的定點FFT演算法

時域抽取法的定點FFT程序3.0 ！

簡介

本示常式序演示了如何計算定點和浮點各個長度和位數的按時間抽取FFT演算法。所有程序在VC6下編譯通過。

最精確的FFT計算最好採用浮點，但浮點對於嵌入式設備計算量太大，有的朋友想採用定點，但網上資料不多，且長度位數要求不一，於是本人製作了這個版本。

此程序可以實現對長度為64，128，256，512，1024，2048的數據進行定點FFT（如果需要其他長度可自行修改）以及浮點版本。

根據處理器處理能力不同，對於數據類型的長度，提供了8位版本，16位版本和32位版本，可用於單片機等不同位數的嵌入式處理器。當然，三個版本主要是給定了不同的處理精度。

此程序採用了運算級間動態縮放捨去最低位，保證了輸出不溢出，且提供了儘可能高的精度。具有自動適應輸入值大小的能力。因此有較高的信噪比。

使用說明

FFT計算內容整合在timefft.h中，可將其加入到圖形環境中或者控制台中，正如示例所示。

有兩個值需要根據環境修改：

長度通過#define N 長度值

決定，長度值可取上文所述六種，主要是考慮嵌入式設備處理能力，一般來講8位機用256點，這樣可以保證下標在8bits以內，16位機可以考慮稍大的點數。

類型數據位數通過：

#define BIT 位數

決定，只能取8，16，32

1.對於圖形界面演示示例，信號源和繪圖程序在WaveAnalysisView.cpp中的void CWaveAnalysisView::OnDraw(CDC* pDC)函數中。

測試數據採用兩個正弦波的疊加，即sin(2*3.14159*i/8)+sin(2*3.14159*i/64);

並將其擴大到最大值之間。

繪圖的版本頻域圖像模值沒有開根號，表示的是模值的平方。

如果時域直流量過大，則會引起頻率的0點值過大，圖像變得縮放的太小，如果為了繪圖方便，可將其歸零。

2.控制台下示例版本是將timefft.h文件獨立出來，和一個控制台結合形成的。 其演示的是一個門函數引發的sinc函數。

定點演算法原理和精度分析

由於定點運算位數為定值，因此不能像浮點一樣進行各種乘法和加法處理，因為這兩種運算會產生溢出。同時，對於旋轉因子的小數模式也需要修整為可以運算的整數模式，這便是FFT定點的兩個需要處理的難點。

本演算法採用了運算量較大但效果較好的動態溢出檢查。也就是說在計算每一級之前，對所有值進行遍歷，如果發現有超過溢出允許輸入值的變數，則對整個數組進行右移縮小處理。這樣就保證了數值較小的輸入不會被過度縮放，而數值很大的數據也能夠通過FFT運算。由於對所有的數組值進行縮放，因此不會造成失真，但是幅度大小會有變化。因此可以理解為定點演算法只求解相對值。在使用過程中，盡量不要疊加過多直流分量，因為這將導致頻率0點處值很大，其它的值被縮放殆盡。影響精度。

下面詳細介紹縮放方法：

一個乘加運算可以理解為：

X1+X2*WN

我們需要找出在什麼情況下輸出會最大，並在設定輸出在此數據類型允許值下時，輸入最大取多少。

由於X1，X2，WN都為複數，且abs(X2)=abs(X2*WN)，因此最壞的情況是，在模一定的情況下，原來X2的實部和虛部的數值經過旋轉因子，被全部轉移到了實部，而虛部為零。或者反之。

而由於當X2的實部和虛部相等時，它的模值最大，假設全部轉移到實部，在加上X1的實部的值就構成了最大的輸出。

我們另這個最大的輸出等於類型最大允許值（8bit時為0x7f，16bit時為0x7fff），由於我們需要對所有輸入做統一的限定值以便程序判斷，所以只能有X1.real=X2.real<max，此值得到的最壞情況下輸出為

Sqrt(X2.real^2+X2.imag^2)+X1.real= Sqrt(2*X2.real^2)+X1.real=(sqrt(2)+1)X1.real

我們將此輸出值限定在類型數據允許的最大整數上，則有

(sqrt(2)+1)X1.real<MAX

其中(sqrt(2)+1)為常數=2.4142

得到X1.real最大不能超過MAX/2.4142，虛部一樣。我們定義此值為OVS

蝶形運算每算完一級，就要進行溢出判斷，如果發現其中任何一個值大於此值，則需要對所有值縮放處理。

這時，我們面臨怎樣縮放的問題，如果採用溢出時所有值/2的方法，那麼，它的輸入端允許的最大信號幅值就是幾乎每級計算完後都需要除以2才通過所有計算。此時

輸入/2*2.4142=輸出

如果只計算一級，那麼輸入不得超過OVS/1.207=52(以8bit為例)，

而每計算一級都要乘以1.207，成為底數大於1等比數列，此時如果計算1024點的FFT，則要保證輸出不溢出，則：

輸入*1.207^log(2,N)=輸入*1.207^10<max

以八位為例max=127（0x7f），則輸入不得超過19，可見，信噪比太低，不能應用。

繼續觀察運算我們發現，最多可能在同一級上溢出兩個bit，我們動態的調整縮放的量級則會改善上述情況。也就是在溢出兩位時除以4，溢出一位時除以2，經過計算機模擬，發現此時只要滿足單極不溢出的條件（如8bit時為52）則總體就不會溢出，但是，付出的代價是兩個bit的精度被處理掉了。因此整體的精度為7bit-2bit（8bit-符號位）=5bit，當然誤差可能由於多次這樣的操作而累積。

程序中，實現時在每次進入新的一級運算之前，將所有值遍歷，如果超出2倍的ovs值則除以4，超出一倍除以2，否則不處理。這樣，乘加運算在數據類型內部不會造成溢出了。

旋轉因子的縮放：由於在計算旋轉因子時，不能用小數，因此，將小數擴大到與位數相當，之後與x[l]相乘，結果再以相同比例縮小，

比如32位版本中：首先將其擴大

wn[k].real=(int)(cos(2*3.14159*k/N)*0x7fffffff);

在計算關於旋轉因子的乘積項時乘完的結果再縮小。

XkWn.real=(((__int64)x[l].real*wn[wi].real)>>31)-(((__int64)x[l].imag*wn[wi].imag)>>31);

XkWn.imag=(((__int64)x[l].imag*wn[wi].real)>>31)+(((__int64)x[l].real*wn[wi].imag)>>31);

(__int64)是由於VC不支持32位乘法時自動切換到64為導致了數據丟失，因此做強制轉換，其它類型不需要。

下載地址：UploadFiles/2009-2/231626252129.rar