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