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