【正文】 為 DIF)法，如果算法是通過逐次分解時間序列 x(n)得到的，這種算法稱為按時間抽取的算法。本文選擇的是按時間抽取的 FFT 算法(DITFFT)。 按時間抽取(DIT)的基2FFT 算法原理設序列 x(m)是長度 N=2m 的有限長序列，按 n 的奇偶把 x(m)分為兩個 N/2點的奇偶子序列：? x(2r ) = x1 (r )?? x(2r + 1) = x2 (r )則可將 DFT 化為：X (k ) = DFT [ x(n)] = ∑ x (n)Wn =0N?12r =0N?12r =0N ?1nkNN ?1r = 0,1,...N ?12N ?1(33)=n =0n偶數(shù)∑ x(n)W nkN+n= 0n奇數(shù)∑ x(n)WnkN 2(2= ∑ x(2r )WN rk + ∑ x(2r + 1)WN r +1) kN?12r =0N?12r =0= ∑ x1 (r )(WN2 )rk + WNk ∑ x2 (r )(WN2 )rk 15 (34)哈爾濱工業(yè)大學工學碩士學位論文2π 2NN )2由于 W = e2N?j=e? j 2π /(= WN / 2 ，故上式可表示為：N?12r =0 rkkkX (k ) = ∑ x1 (r )WN / 2 + WN ∑ x2 (r )WNrk/2 = X 1 (k ) + WN X 2 (k )r =0N?12(35)式中 X1(k)和 X2(k)分別是 x1(r)和 x2(r)的 N/2 點 DFT。由此我們可以看到，一個 N 點的 DFT 已分解為兩個 N/2 點的 DFT。這兩個N/2 點的 DFT 再按照式(35)組合成一個 N 點的 DFT。應該看到，X1(k)和 X2(k)只有 N/2 個點，即 k=0,1,…,N/21，而 X(k)卻有 N 個點，即 k=0,1,…,N1，故式得到的只是 X(k)的前一半的結果，要用 X1(k)和 X2(k)來表達全部的 X(k)值，還必須利用系數(shù)的周期性：W這樣可得到： X 1 (同理可得： X 2 (kNN?12rkN /2r (k +N )2N?12r =0= WN /2(36)N+ k ) = ∑ x1 (r )W2r =0 Nr (k + ) 2N /2 rk= ∑ x1 (r )WN / 2 = X 1 (k )N + k ) = X 2 (k )2(N+k )2 kk= WNN / 2WN = ?WN ，這樣可得到 X (k ) 后半部分再考慮到 W 的性質： WN的表達式。綜上， X (k ) 表達為前后兩部分為： k? X (k ) = X 1 (k ) + WN X 2 (k )??N X (k + ) = X 1 (k ) ? WNk X 2 (k )??2k = 0,1,...,N ?12(37)可見，只要求出 0~(N/2－1)區(qū)間內的所有 X1(k)和 X2(k)值，就可求出 0~(N－1)區(qū)間內的所有 X(k)值，顯然節(jié)省了運算量。式(37)的運算可以用下圖的蝶形信號流圖(又稱蝶形運算單元)來表示。X 1 (k )X 1 ( k ) + WNk X 2 ( k )X 2 (k ) kWNX 1 ( k ) ? WNk X 2 (k )圖 31 時間抽取法蝶形運算流圖 16 哈爾濱工業(yè)大學工學碩士學位論文 按時間抽取的 FFT 算法的特點原位運算（同址運算）原位運算可節(jié)省存儲單元，降低設備成本。蝶形運算是很有規(guī)律的，其每級（每列）計算都是有 N/2 個蝶形運算構成，每一蝶形結構完成下述迭代運算： r? X m (k ) = X m ?1 (k ) + WN X m?1 ( j )?r? X m ( j ) = X m ?1 (k ) ? WN X m?1 ( j )(38)式(38)中，m 表示第 m 列迭代，k，j 為數(shù)據所在行數(shù)。式(38)的蝶形運算用圖 31 來表示，由一次復數(shù)乘法和兩次復數(shù)加（減）法組成。可以看出，某一列的任何兩個節(jié)點 k 和 j 的節(jié)點變量進行蝶形運算后，得到結果為下一列 k和 j 兩個節(jié)點的節(jié)點變量，和其他節(jié)點變量無關，因而可采用原位運算，也就是蝶形的兩個輸出值仍放回蝶形的兩個輸入所在的存儲器中。倒位序規(guī)律當運算完成后，F(xiàn)FT 的輸出 X(k)按正常順序排列在存儲單元中，但這時輸入 x(n)卻不是按自然順序存儲的，這正是由于對 x(n)做奇、偶分開所產生的。對于 N=8，設序列的自然順序號 I 用二進制表示為：(n2 n1 n0)，則其倒位序J 對應的二進制數(shù)就是：(n0 n1 n2)。自然順序數(shù) I 增加 1，是在順序數(shù)的二進制最低位加 1，向左進位，而倒序數(shù) J 則是在二進制數(shù)最高位加 1，向右進位。用這種算法，可以由當前任一倒序值求得下一個倒序值。N=8 時把按自然順序存放在存儲器單元中的數(shù)據，換成 FFT 原位運算流圖所要求的倒位序變址功能如圖 32 所示，當 I=J 時，不必調換；當 I≠J 時，需要將原來存放數(shù)據 x(I)的存儲單元內調入數(shù)據 x(J)，而存放 x(J)的存儲單元內調入x(I)。為了避免把已掉換過的數(shù)據再次調換，保證只調換一次（否則又回到原狀），我們只需看 J 是否不 I 小。若 JI，意味著此 x(I)在前面已和 x(J)調換過，不必再調換了； JI，若才將兩存儲單元內容互換。這樣就得到輸入所需的倒位序列。圖 32 N=8 倒位序的變址處理 17 哈爾濱工業(yè)大學工學碩士學位論文蝶形運算兩節(jié)點的“距離”以圖 32 的 8 點的 FFT 為例，其輸入是倒序的，輸出是自然順序的。其第一級（第一列）每個蝶形運算的兩節(jié)點間“距離”為 1，第二級每個蝶形運算的兩節(jié)點間“距離”為 2，第三極每個蝶形運算的兩節(jié)點間“距離”為 4，依次類推，對 N=2M 點 FFT，當輸入為倒位序，輸出為正常順序時，第 m 級運算每個蝶形的兩節(jié)點“距離”為 2m1。 r WN 因子的確定對于第 m 級運算，一個 DFT 蝶形運算的兩節(jié)點“距離” 2m1，為因而式(38)可以寫成： r? X m (k ) = X m ?1 (k ) + X m ?1 (k + 2 m?1 )WN?m ?1m ?1r? X m (k + 2 ) = X m ?1 (k ) ? X m?1 (k + 2 )WN r為了完成式(39)的運算，還必須知道系數(shù) WN 的變換規(guī)律：(39)(1) 把上式中蝶形運算兩節(jié)點中的第一個節(jié)點標號值 k，表示成 M 位二進制數(shù)（N=2M）；(2) 把此二進制數(shù)乘上 2Mm，即將此 M 位二進制數(shù)左移 Mm 位，（m 表示第 m 級運算），右邊空出來的位置補零，此數(shù)即為所求 r 的二進制數(shù)。由此可得出， W 因子最后一列有 N/2 種，順序為 W ， W ，…， WNrN0N1N(N?1)2，其余可類推。 實序列 FFT 的高效算法在發(fā)現(xiàn)了 FFT 結果的線性和奇偶對稱性之后，人們認識到兩個 N 點實序列的 FFT 可以使用一個 N 點的復序列的 FFT 來實現(xiàn)，并提出了使用 N 點復序列的 FFT 求取 2N 點的實序列 FFT 的方法。由于實際工作中，A/D 采樣的值一般都是實數(shù)，即輸入數(shù)據 x(n)為實序列，而一般的 FFT 算法都是基于復數(shù)，即每個輸入數(shù)都含有實部和虛部。如不采取特殊措施，往往是把 x(n)視為一個虛部為零的復序列，這樣會增加運算量。我們可以利用 FFT 的性質將 2N 個點的實輸入序列組合成一個 N 點的復序列，然后再對復序列進行 N 點的 FFT 運算，最后再由 N 點的復數(shù)輸出拆成 2N 點的復數(shù)序列， 2N 點的復數(shù)序列與原始 2N 點的實數(shù)序列的 DFT 輸出一致。這使用這種方法，在組合輸入和拆散輸出的操作中，F(xiàn)FT 運算量減半，這樣利用實數(shù) FFT 18 哈爾濱工業(yè)大學工學碩士學位論文算法來計算實輸入序列的 DFT 的速度幾乎是一般復 FFT 算法的兩倍。下面就詳細的分析這種針對實輸入數(shù)據的高效 FFT 算法，其基本思路是將2N 點的實序列分成奇偶兩個序列，用 N 點復數(shù) FFT 運算代替 2N 點的實數(shù) FFT運算。對于一個 2N 點實序列 a(n)，其中 0≤n≤2N1，要計算的 2N 點傅里葉變換Xa(m)。把 a(n)的奇偶序列值恰當?shù)卮嫒?N 點復序列 x(n)，即x(0)=a(0)+ja(1)x(1)=a(2)+ja(3)x(2)=a(4)+ja(5) Mx(N－1)=a(2N－2)+ja(2N－1)(310)把式(310)中的復序列用于 N 點復序列 FFT 算法，就可以得到 FFT 結果X(m)=Xr(m)+jXi(m)，其中 m 從 0 到 N－1。為了從 X(m)中提取出想要的 2N 點實序列 FFT 結果 Xa(m)=Xa,real(m)+jXa,imag(m)，定義以下關系式：X r+ (m) = X r ( m) + X r ( N ? m) 2 X ( m) ? X r ( N ? m)X r? (m) = r 2 X (m) + X i ( N ? m)X i+ (m) = i 2 X ( m) ? X i ( N ? m)X i? (m) = i 2(311)(312)(313)(314)把式(311)~ (314)的值代入下面的式子就可以求出最后的結果 Xa(m)的實部和虛部： πmπm ) ? X i+ (m) ? sin() ? X r? (m) NN πmπmX a ,imag (m) = X i? (m) ? sin() ? X i+ (m) ? cos() ? X r? (m) NNX a , real (m) = X r+ (m) + cos((315)(316)這里需要注意，原序列 a(n)中的 n 變化范圍是 0 到 2N－1，而 N 點 FFT 結果中 m 的變化范圍是 0 到 N－1。這個算法只能得到 2N 點實序列 FFT 結果的前半部分 Xa(0)~Xa(N－1)。因為此算法要求輸入序列 a(n)為實數(shù)，所以 Xa(N)~Xa(2N－1)是 Xa(0)~Xa(N－1)的復共軛，不需要重新計算。圖 33 給出了 2N 點實序列 19 哈爾濱工業(yè)大學工學碩士學位論文FFT 算法的計算流程圖。X (m) X (m)X i+ (m) X i? (m)+r?rX a( m ) =X a ,real ( m ) + jX a,imag ( m )X a (m)圖 33 2N 點實序列 FFT 算法流程這種改進算法可以使一般的 FFT 運算量減少差不多一半，表 31 給出了不同長度的實數(shù)序列采用 DFT、FFT、改進 FFT 算法的運算量（實數(shù)乘法和實數(shù)加法），直觀的說明了改進算法帶來的好處。表 31 不同算法運算量的比較序列長度 128 256 5121024 DFT 算法 乘法加法 327683264013107213081652428852377620971522096128 FFT 算法 乘法加法 17922688 40966144 9216138242048030720 改進 FFT 算法 乘法加法 10241664 23043712 512081921126417920 諧波分析中 FFT 算法存在的問題及解決辦法在應用 DFT 分析電力信號的實際問題中會遇到下列幾個問題：混疊現(xiàn)象、頻譜泄露、柵欄效應[16]。這些問題對系統(tǒng)的影響是不同的，本節(jié)將對以上問題作逐一分析，并提出本系統(tǒng)中的解決方法?；殳B現(xiàn)象根據乃奎斯特抽樣定理：一個頻譜有限的信號 f(t)，如果頻譜只占據－Ωk~－Ωk 的范圍，則信號 f(t)可以用等間隔的采樣值來表示，而抽樣頻率必須大于信號的最高頻率的兩倍即 2Ωk。然而對大多數(shù)信號而言，頻