數字訊號處理公式變程式（一）——DFT、FFT
之前搞了一些數字訊號處理演算法程式設計（OC），一直沒來得及整理，現在整理一下。陸續會更新，包括FFT、巴特沃斯濾波器（高低帶通、高低帶阻）、資料差值（線性、sinc、三次樣條*）、資料壓縮（等距、平均、峰值檢測）和模仿matlab的STFT功能（spectrogram函式三維繪圖）。
注：我比較喜歡使用matlab（也喜歡自己修改matlab的原始碼做測試，所以重灌了好幾次了，囧。。。），有部分參考matlab的演算法進行（或者直接翻譯為OC），演算法的執行結果經常會跟matlab運算結果進行比較，演算法只做學習用，轉載請註明。
另外可能會有不足或者理解偏差的地方，路過的高人請不吝賜教。
好啦，進入正題。
---------------------------------------------------------------------------------------
在數字世界中FFT的意義不言而喻（我曾轉載一篇文章有提到：http://blog.csdn.net/shengzhadon/article/details/40539101），這裡就不再贅述了。FFT（快速傅立葉變換）是DFT的一種特殊情況，就是當運算點的個數是2的整數次冪的時候進行的運算（不夠用0補齊），那就先從DFT開始吧。
一、DFT（本部分就是翻譯公式）
定義（來自百科）：離散傅立葉變換（Discrete Fourier Transform，縮寫為DFT），是傅立葉變換在時域和頻域上都呈離散的形式，將訊號的時域取樣變換為其DTFT的頻域取樣。在形式上，變換兩端（時域和頻域上）的序列是有限長的，而實際上這兩組序列都應當被認為是離散週期訊號的主值序列。即使對有限長的離散訊號作DFT，也應當將其看作其週期延拓的變換。在實際應用中通常採用快速傅立葉變換計算DFT。
DFT的公式是：設x(n)為M點的有限長序列，即在0≤n≤M-1內有值，則定義x(n)的N點（N≥M。當N＞M時，補N-M個零值點），離散傅立葉變化定義為
其中，
①週期性：W(N, nk) = W(N, (n+rN)k) = W(N, n(k+rN))，其中r問整數。
②共軛對稱性：(W(N, nk))* = W(N, -nk)。
③可約性：W(N, nk) = W(mN, mnk)，W(N, nk) = W(N/m, nk/m)，其中m為整數，N/m為整數。
④特殊值：
W(N, N/2) = -1； W(N, (k+N/2)) = -W(N, k)； W(N, (N-k)n) = W(N, (N-n)k) = W(N, -nk)。
所以，在計算旋轉因子的過程中可以適當的使用特殊值來提高運算的效率。
在計算DFT時，如果資料的點數夠計算的點數，則擷取計算點數長的資料進行DFT運算，否則將資料點個數補0至計算點個數，然後進行計算，舉例如下（舉例只是為了說明問題，沒有按照程式語言的書寫格式）。
比如：資料點陣列為 dataArray = {1， 2， 3， 4，5}，可以看出資料長度為5。
如果要求做4點DFT運算，則只需擷取前4個數作為運算陣列進行運算即可，即為{1， 2， 3， 4}；
如果要求做8點DFT運算，則需在原陣列後補三個0，使長度為8後再進行計算，計算陣列為{1， 2， 3， 4，5，0，0，0}。
旋轉因子計算程式碼如下：
注：①傳入值格式為W(N, p)；②為了提高效率，當旋轉因子p = 0時，直接返回 result = 1+0*j = 1；③返回值為負數形式；④利用尤拉公式求負數的指數運算。
/*======================================================================
* 方法名： twiddleFactor
* 方法功能：求FFT計算過程中的旋轉因子 - W(N,p)
* 說明： Euler（尤拉）公式 exp(iθ) = cos(θ) - i*sin(θ)
*
* 變數名稱：
* N - FFT計算點數
* p - 旋轉因子的階數
*
* 返回值： 旋轉因子的複數表示
*=====================================================================*/
+ (ComplexStruct *)twiddleFactor:(int)N andP:(int)p
{
ComplexStruct *twiddle = [[ComplexStruct alloc] init];
if(p==0)
{
[twiddle setReal:1 andSetImg:0];
}
else
{
float tempx = PI_x_2 * p/N;
float real = cos(tempx); //計算複數的實部與虛部
float img = -1 * sin(tempx);
[twiddle setReal:real andSetImg:img]; //呼叫賦值方法
}
return twiddle;
}DFT流程圖和程式碼如下：
/*======================================================================
* 方法名： dft
* 方法功能：計算陣列的DFT，非2的整數次冪的FFT
* N-1
* 公式說明：X(k) = ∑ [x(n)*W(N,nk)]，其中，W(N,nk)為旋轉因子，k = 0,1,...,N-1
* n=0
*
*
* yVector - 原始資料
* length - 原始資料長度
* N - FFT計算點數
* fftYreal - FFT後的實部
* fftYImg - FFT後的虛部
*
* 返回值： 是否成功的標誌，若成功返回true，否則返回false
*=====================================================================*/
+ (BOOL)dft:(float *)yVector andOriginalLength:(NSInteger)length andFFTCount:(NSInteger)N andFFTReal:(float *)fftYReal andFFTYImg:(float *)fftYImg
{
// 標誌
BOOL isFFTOK = false;
// 旋轉因子指數
NSInteger p;
// N點運算的原始資料
float yVectorN[N];
for (int i = 0; i < N; i++)
{
yVectorN[i] = 0.0;
}
// 定義計算過程中用到的複數變數
ComplexStruct *tempFFTY, *yFFT1, *yFFT2;
// 初始化複數變數
yFFT1 = [[ComplexStruct alloc]init];
yFFT2 = nil;
// 判斷點數是否夠FFT運算點數，並確定最終用於計算的點的值
if (length <= N)
{
// 如果N至少為length，先把yVector全部賦值
for (NSInteger i = 0; i < length; i++)
{
yVectorN[i] = yVector[i];
}
// 如果 N > length，後面補零
if (length < N)
{
for (NSInteger i = length; i < N; i++)
{
yVectorN[i] = 0.0;
}
}
}
else
{
// 如果 N < length 擷取相應長度的資料進行運算
for (NSInteger i = 0; i < N; i++)
{
yVectorN[i] = yVector[i];
}
}
// 計算結果初始化
for (NSInteger i = 0; i < N; i++)
{
fftYReal[i] = 0.0;
fftYImg[i] = 0.0;
}
// 運用公式計算
for (NSInteger i = 0; i < N; i++)
{
for (NSInteger j = 0; j < N; j++)
{
p = i*j;
tempFFTY = [FFT twiddleFactor:N andP:p];
[yFFT1 setReal:yVectorN[j] andSetImg:0.0];
yFFT2 = [ComplexStruct complexMul:yFFT1 and:tempFFTY];
fftYReal[i] += [yFFT2 getReal];
fftYImg[i] += [yFFT2 getImg];
}
}
return isFFTOK;
}二、FFT
1.DFT的運算量
複數乘法次數為N*N，複數加法次數為N(N-1)，若N>>1，則這兩者都近似為N*N，它隨N增大為急速增大。
改進途徑：①利用旋轉因子性質減小計算量；②由於運算量和N*N成正比，因而可將N點DFT分解成小點數的DFT，以減少運算量（點數越小，計算量越小）；③改為用FFT計算，複數乘法次數為(N/2)*log(2,N)。
2.FFT可分為按照時間抽選的基-2演算法（庫利-圖基演算法DIT-FFT）和按頻率抽選的基-2演算法（桑德-圖基演算法DIF-FFT）。本文采用DIT-FFT演算法。
3.FFT計算原理及流程圖
FFT的計算要求點數必須為2的整數次冪，如果點數不夠用0補齊。例如計算{2，3，5，8，4}的8點FFT，需要補3個0後進行計算，如果計算該陣列的5點FFT，則先計算8點FFT後擷取前5個值即可（不提倡）。
（1）原理
公式推導
設序列x(n)點數為N=2^L，L為整數。將N=2^L，先按照n的奇偶分為兩組，其中r = 0, 1, ..., N/2-1
x(2r) = x1(r)
x(2r-1) = x2(r)
則可將DFT化為（式1）
由上式可以看出一個N點的DFT可以分為兩個N/2點的DFT，按照上式右組合成N點DFT。但是這裡的x1(r)、x2(r)以及X1(k)、X2(k)都是N/2點的序列，即r,k滿足r,k=0, 1, ..., N/2-1。而X(k)卻有N點，上式計算的只是X(k)的前半項數的結果，因此還需要計算後半項的值。
將①②③代入（式1），就可將X(k)表達為前後兩部分：
前半部分，X(k)，當k=0, 1, ..., N/2-1
後半部分，X(k)，當k=N/2, ..., N-1
蝶形運算說明
蝶形運算，符號表示如下圖所示：
所以，FFT的蝶形運算圖表示為（8點，2^3 = 8，運算級數最大為L=3）
在蝶形運算中變化規律由W(N, p)推導，其中N為FFT計算點數，J為下角標的值
L = 1時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0;
L = 2時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1;
L = 3時，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3;
所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1
又因為2^L = 2^M*2^(L-M) = N*2^(L-M)，這裡N為2的整數次冪，即N=2^M，
W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))
所以，p = J*2^(M-L)，此處J = 0, 1, ..., 2^(L-1)-1，當J遍歷結束但計算點數不夠N時，J=J+2^L，後繼續遍歷，直到計算點數為N時不再迴圈。
舉例：N=8點的FFT計算
當L=2時，J = 0, 1兩個值，因此p = J*2^(M-L) = 0, 2兩個值，即旋轉因子有兩個值W(8, 0)和W(8, 2)，計算中兩行之間的距離B = 2^(L-1)=2^(2-1)=2，
代入J=0, 1可求得X(0)、X(0+2)和X(1)、X(1+2)，即可求出第2級蝶形運算的X(0), X(1), X(2), X(3)，也就是求出一半，此時J加步進2^L=4，即J=J+2^L=4, 5，
再代入J=4, 5可求出X(4), X(4+2)和X(5), X(5+2)，即可求出第2級蝶形運算的X(4), X(5), X(6), X(7)，已經全部求出，J迴圈結束。
二進位制倒序說明
前面數的排列順序是進行二進位制倒序後的排序。二進位制倒序是指將某數轉化為二進位制表示，將最高位看做最低位、次高位看做次低位...以此類推，計算後的紙進行排序。例如3的二進位制表示為011b（3=0*2^2+1*2+1），二進位制倒序值為6（011b，最高位看做最低位...即6=0+1*2+1*2^2）
即0到7的二進位制排序是：
0, 4, 2, 6, 1, 5, 3, 7
000→000 001→100 010→010 011→110 100→001 101→101 110→011 111→111
（2）流程圖及程式碼
①FFT運算總流程圖，包括
②二進位制倒序程式碼（過程見程式碼註釋，可用其他方式實現）
// 二進位制倒序說明：
// 二進位制倒序是指將數二進位制表示後，把最高位當做最低位（以此類推）排序，最高位向次高位進位
// 如倒序後 0 = 000，1 = 100（最高位+1），2 = 010（最高位向次高位進1）
//------------------------------------------------------------------------
// 倒序前 高位->低位 低位->高位 倒序後
//
// 0 0 0 0 --- 0 0 0 0
// 1 0 0 1 --- 1 0 0 4
// 2 轉化為二進位制 0 1 0 --- 0 1 0 轉化為十進位制 2
// 3 ===========> 0 1 1 --- 1 1 0 ===========> 6
// 4 1 0 0 --- 0 0 1 1
// 5 1 0 1 --- 1 0 1 5
// 6 1 1 0 --- 0 1 1 3
// 7 1 1 1 --- 1 1 1 7
//------------------------------------------------------------------------
/**
* @brief 對陣列vector進行二進位制倒序處理，按位計算
* 實現原理：將原陣列下角標倒序後賦值給倒序後陣列
* @param[in] vector 原陣列
* @param[in] N 陣列的長度
* @param[out] inverseOrderVector 倒序後的陣列
*/
+ (void)inverseOrder:(float *)vector andN:(NSInteger)N andInverseOrderVector:(float *)inverseOrderVector
{
// 求對數，即求出二進位制表示時N的位數
NSInteger k = log2(N);
// 按位處理
for (NSInteger i = 0; i < N; i++)
{
// 下角標號
NSInteger temp_i = i;
// 計算後的下角標
NSInteger foot = 0;
// 計算過程中，當前所在的位
NSInteger j = 0;
// 第j位上的數
NSInteger temp_foot = 0;
// 每次取低位，變為對應高位，迴圈求和
while (temp_i)
{
// 取出最後一位
temp_foot = temp_i % 2;
// 計算將最後一位轉化為高位的十進位制值，並與之前值累加
foot = foot + temp_foot * powf(2, (k - 1 - j));
// 將數右移一位，最右位捨棄
temp_i = temp_i * 0.5;
// 下一位
j++;
}
// 根據下角標取值
inverseOrderVector[i] = vector[foot];
}
}③根據（1）中推導的內容，FFT的流程圖可化為
說明：蝶形運算中又三層迴圈
第一層（最外層），完成M次迭代過程，即算出A0(k), A1(k), ..., Am(k)，其中k=0, 1, ..., N；A2(k)為蝶形運算第2級的結果，如A0(k)=x(k), Am(k)=X(k)；
第二層（中間層），完成旋轉因子的變化，步進為2^L；
第三層（最裡層），完成相同旋轉因子的蝶形運算。
FFT實現的程式碼如下：
/*======================================================================
* 方法名： fft
* 方法功能：計算陣列的FFT，運用蝶形運算
*
* 變數名稱：
* yVector - 原始資料
* length - 原始資料長度
* N - FFT計算點數
* fftYreal - FFT後的實部
* fftYImg - FFT後的虛部
*
* 返回值： 是否成功的標誌，若成功返回true，否則返回false
*=====================================================================*/
+ (BOOL)fft:(float *)yVector andOriginalLength:(NSInteger)length andFFTCount:(NSInteger)N andFFTReal:(float *)fftYReal andFFTYImg:(float *)fftYImg
{
// 確保計算時時2的整數冪點數計算
NSInteger N1 = [self nextNumOfPow2:N];
// 定義FFT運算是否成功的標誌
BOOL isFFTOK = false;
// 判斷計算點數是否為2的整數次冪
if (N != N1)
{
// 不是2的整數次冪，直接計算DFT
isFFTOK = [self dft:yVector andOriginalLength:length andFFTCount:N andFFTReal:fftYReal andFFTYImg:fftYImg];
// 返回成功標誌
return isFFTOK;
}
// 如果計算點數位2的整數次冪，用FFT計算，如下
// 定義變數
float yVectorN[N1]; // N點運算的原始資料
NSInteger powOfN = log2(N1); // N = 2^powOfN，用於標記最大運算級數(公式中表示為：M)
NSInteger level = 1; // 運算級數（第幾次運算），最大為powOfN，初值為第一級運算(公式中表示為：L)
NSInteger lineNum; // 行號，倒序排列後的蝶形運算行號(公式中表示為：k)
float inverseOrderY[N1]; // yVector倒序x
NSInteger distanceLine = 1; // 行間距，第level級運算每個蝶形的兩個節點距離為distanceLine = 2^(L-1)(公式中表示為：B)
NSInteger p; // 旋轉因子的階數，旋轉因子表示為 W(N, p)，p = J*2^(M-L)
NSInteger J; // 旋轉因子的階數，旋轉因子表示為 W(2^L, J)，J = 0, 1, 2,..., 2^(L-1) - 1 = distanceLine - 1
float realTemp, imgTemp, twiddleReal, twiddleImg, twiddleTheta, twiddleTemp = PI_x_2/N1;
NSInteger N_4 = N1/4;
// 判斷點數是否夠FFT運算點數
if (length <= N1)
{
// 如果N至少為length，先把yVector全部賦值
for (NSInteger i = 0; i < length; i++)
{
yVectorN[i] = yVector[i];
}
if (length < N1)
{
// 如果 N > length 後面補零
for (NSInteger i = length; i < N1; i++)
{
yVectorN[i] = 0.0;
}
}
}
else
{
// 如果 N < length 擷取相應長度的資料進行運算
for (NSInteger i = 0; i < N1; i++)
{
yVectorN[i] = yVector[i];
}
}
// 呼叫倒序方法
[self inverseOrder:yVectorN andN:N1 andInverseOrderVector:inverseOrderY];
// 初始值
for (NSInteger i = 0; i < N1; i++)
{
fftYReal[i] = inverseOrderY[i];
fftYImg[i] = 0.0;
}
// 三層迴圈
// 第三層（最裡）：完成相同旋轉因子的蝶形運算
// 第二層（中間）：完成旋轉因子的變化，步進為2^level
// 第一層（最外）：完成M次迭代過程，即計算出x(k) = A0(k), A1(k),...,Am(k) = X(k)
// 第一層迴圈
while (level <= powOfN)
{
distanceLine = powf(2, level - 1); // 初始條件 distanceLine = 2^(level-1)
J = 0;
NSInteger pow2_Level = distanceLine * 2; // 2^level
NSInteger pow2_NSubL = N1/pow2_Level; // 2^(M-L)
// 第二層迴圈
while (J < distanceLine)
{
p = J * pow2_NSubL;
lineNum = J;
NSInteger stepCount = 0; // J運算的步進計數
// 求旋轉因子
if (p==0)
{
twiddleReal = 1.0;
twiddleImg = 0.0;
}
else if (p == N_4)
{
twiddleReal = 0.0;
twiddleImg = -1.0;
}
else
{
// 計算尤拉公式中的θ
twiddleTheta = twiddleTemp * p;
// 計算複數的實部與虛部
twiddleReal = cos(twiddleTheta);
twiddleImg = -1 * sin(twiddleTheta);
}
// 第三層迴圈
while (lineNum < N1)
{
// 計算下角標
NSInteger footNum = lineNum + distanceLine;
/*---------------------------------------
* 用複數運算計算每級中各行的蝶形運算結果
* X(k) = X(k) + X(k+B)*W(N,p)
* X(k+B) = X(k) - X(k+B)*W(N,p)
*---------------------------------------*/
realTemp = fftYReal[footNum] * twiddleReal - fftYImg[footNum] * twiddleImg;
imgTemp = fftYReal[footNum] * twiddleImg + fftYImg[footNum] * twiddleReal;
// 將計算後的實部和虛部分別存放在返回陣列中
fftYReal[footNum] = fftYReal[lineNum] - realTemp;
fftYImg[footNum] = fftYImg[lineNum] - imgTemp;
fftYReal[lineNum] = fftYReal[lineNum] + realTemp;
fftYImg[lineNum] = fftYImg[lineNum] + imgTemp;
stepCount += pow2_Level;
// 行號改變
lineNum = J + stepCount;
}
// 旋轉因子的階數變換，達到旋轉因子改變的效果
J++;
}
// 運算級數加一
level++;
}
isFFTOK = true;
return isFFTOK;
}
三、DFT、FFT測試結果
本文開頭就已經說了，我個人比較喜歡使用matlab，所以我程式的對比物件當然是matlab。給定t1[8] = {0,1,2,3,4,5,6,7};signal=sin(t1);對signal進行8、16、4、13、6點FFT/DFT計算輸出的結果和給定128點資料進行128點FFT計算結果如下圖所示：
複製一遍：
FFT，8點計算結果
0.553733 + 0.000000i 2.394647 - 2.097012i -1.386684 + 0.915560i -0.881042 + 0.280414i -0.807574 + 0.000000i -0.881042 - 0.280414i -1.386684 - 0.915560i 2.394647 + 2.097012i
FFT，16點計算結果
0.553733 + 0.000000i 1.431957 + 0.493527i 2.394647 - 2.097012i -1.786256 - 2.899550i -1.386684 + 0.915560i 0.105162 - 0.495158i -0.881042 + 0.280414i 0.249137 - 0.129291i -0.807574 + 0.000000i 0.249137 + 0.129291i -0.881042 - 0.280414i 0.105162 + 0.495158i -1.386684 - 0.915560i -1.786256 + 2.899551i 2.394647 + 2.097012i 1.431957 - 0.493527i
FFT，4點計算結果
1.891888 + 0.000000i -0.909297 - 0.700351i -0.073294 + 0.000000i -0.909297 + 0.700351i
FFT，13點計算結果
0.553733 + 0.000000i 1.898167 + 0.288124i 0.803761 - 3.465451i -2.328951 + 0.137430i 0.200269 - 0.367077i -0.688243 + 0.477332i -0.161869 - 0.528187i -0.161868 + 0.528186i -0.688245 - 0.477332i 0.200271 + 0.367075i -2.328952 - 0.137428i 0.803762 + 3.465451i 1.898167 - 0.288123i
FFT，6點計算結果
0.176162 + 0.000000i -0.276095 - 3.002073i 0.123599 - 0.116304i 0.128828 + 0.000000i 0.123600 + 0.116302i -0.276092 + 3.002073i
FFT，128點計算結果
0.699691 + 0.000000i 0.704772 + 0.010352i 0.719298 + 0.021364i 0.739474 + 0.033697i 0.749607 + 0.047746i 0.679161 + 0.061398i 0.105405 + 0.050647i -6.945906 - 0.423135i 54.641422 + 5.062080i -57.156776 - 6.630915i 9.187446 + 1.338868i 0.342783 + 0.093026i -0.288066 - 0.020889i -0.353394 - 0.043726i -0.328844 - 0.047995i -0.290196 - 0.047167i -0.253328 - 0.044825i -0.221405 - 0.042109i -0.194470 - 0.039423i -0.171855 - 0.036897i -0.152811 - 0.034570i -0.136679 - 0.032440i -0.122920 - 0.030495i -0.111105 - 0.028715i -0.100891 - 0.027084i -0.092007 - 0.025581i -0.084229 - 0.024193i -0.077390 - 0.022905i -0.071342 - 0.021706i -0.065971 - 0.020587i -0.061182 - 0.019537i -0.056895 - 0.018550i -0.053044 - 0.017619i -0.049575 - 0.016737i -0.046441 - 0.015901i -0.043601 - 0.015106i -0.041024 - 0.014347i -0.038679 - 0.013621i -0.036542 - 0.012927i -0.034593 - 0.012260i -0.032808 - 0.011617i -0.031181 - 0.010997i -0.029686 - 0.010398i -0.028318 - 0.009818i -0.027066 - 0.009256i -0.025917 - 0.008710i -0.024865 - 0.008177i -0.023903 - 0.007660i -0.023024 - 0.007153i -0.022220 - 0.006657i -0.021489 - 0.006172i -0.020825 - 0.005695i -0.020224 - 0.005227i -0.019682 - 0.004767i -0.019197 - 0.004313i -0.018766 - 0.003867i -0.018383 - 0.003419i -0.018055 - 0.002986i -0.017772 - 0.002551i -0.017534 - 0.002121i -0.017342 - 0.001693i -0.017193 - 0.001268i -0.017087 - 0.000845i -0.017025 - 0.000423i -0.017003 + 0.000000i -0.017025 + 0.000422i -0.017087 + 0.000845i -0.017193 + 0.001268i -0.017342 + 0.001693i -0.017534 + 0.002121i -0.017772 + 0.002551i -0.018055 + 0.002986i -0.018383 + 0.003425i -0.018766 + 0.003866i -0.019197 + 0.004314i -0.019682 + 0.004767i -0.020224 + 0.005227i -0.020825 + 0.005695i -0.021489 + 0.006172i -0.022220 + 0.006656i -0.023024 + 0.007153i -0.023902 + 0.007659i -0.024865 + 0.008177i -0.025917 + 0.008710i -0.027066 + 0.009256i -0.028318 + 0.009818i -0.029686 + 0.010398i -0.031179 + 0.010998i -0.032808 + 0.011617i -0.034592 + 0.012259i -0.036542 + 0.012927i -0.038679 + 0.013621i -0.041024 + 0.014347i -0.043601 + 0.015106i -0.046441 + 0.015901i -0.049574 + 0.016737i -0.053044 + 0.017619i -0.056894 + 0.018549i -0.061182 + 0.019537i -0.065971 + 0.020587i -0.071342 + 0.021706i -0.077390 + 0.022905i -0.084229 + 0.024193i -0.092005 + 0.025579i -0.100891 + 0.027084i -0.111107 + 0.028715i -0.122921 + 0.030496i -0.136679 + 0.032441i -0.152811 + 0.034570i -0.171855 + 0.036897i -0.194470 + 0.039423i -0.221401 + 0.042109i -0.253328 + 0.044825i -0.290195 + 0.047166i -0.328844 + 0.047996i -0.353394 + 0.043726i -0.288066 + 0.020889i 0.342783 - 0.093026i 9.187446 - 1.338871i -57.156776 + 6.630923i 54.641422 - 5.062086i -6.945903 + 0.423136i 0.105404 - 0.050648i 0.679161 - 0.061398i 0.749607 - 0.047746i 0.739474 - 0.033698i 0.719298 - 0.021364i 0.704768 - 0.010354i
測試結果顯示：iOS程式與MATLAB程式運算結果一致，但是MATLAB運算效率高。以下是多次比較求平均後的時間值
為了更直觀的觀察，將其他的訊號（幾個訊號疊加後的訊號）通過fft計算後並作頻域處理匯入matlab，繪製圖如下：
===================================================================
OK，FFT隨筆一份，方便自己以後查閱。其中一定會有很多不足，望路過的大神訂正。謝過