數(shù)字信號(hào)處理CH9 DFT計(jì)算實(shí)現(xiàn)方法-FFT

CH9DFT計(jì)算實(shí)現(xiàn)方法—FFT9.1DFT運(yùn)算量估計(jì)及算法分析9.2Goertzel算法9.3按時(shí)間抽取的FFT算法9.4按頻率抽取的FFT算法9.5FFT實(shí)現(xiàn)問題討論9.6有限寄存器長(zhǎng)度的影響19.1DFT運(yùn)算量估計(jì)及算法分析算法優(yōu)劣的衡量與應(yīng)用相關(guān)，工程中經(jīng)?？紤]的有：精度、信噪比實(shí)時(shí)性、收斂性存儲(chǔ)量、硬件成本運(yùn)算量：乘法與加法次數(shù)各項(xiàng)要求有時(shí)矛盾，研究人員在實(shí)現(xiàn)中要做好和諧工作2例：三角函數(shù)的實(shí)現(xiàn)實(shí)時(shí)性要求很高時(shí)用查表法，如雷達(dá)PPI極坐標(biāo)到光柵掃描直角坐標(biāo)的轉(zhuǎn)換通用性要求較好時(shí)常用冪級(jí)數(shù)展開法，如一般計(jì)算機(jī)程序的函數(shù)庫(kù)輸出有規(guī)則時(shí)可用遞推法，如波形發(fā)生器運(yùn)算量實(shí)時(shí)性存儲(chǔ)量靈活性查表法小極好很大一般Z變換遞推法中很好小差冪級(jí)數(shù)展開法大較好較小好3DFT運(yùn)算量估計(jì)X[k]=x[n]WNkn={Re[x[n]]Re[WNkn]-Im[x[n]]Im[WNkn]+jIm[x[n]]Re[WNkn]+jRe[x[n]]Im[WNkn]}k=0,1,…N-1復(fù)數(shù)N點(diǎn)DFT運(yùn)算：N2次復(fù)乘，N(N-1)次復(fù)加實(shí)數(shù)N點(diǎn)DFT運(yùn)算：4N2次實(shí)乘，N(4N-2)次實(shí)加4DFT算法分析算法分析的目的是減少乘法與加法的次數(shù)Re[WNkn]=cos(2πkn/N)，Im[WNkn]=-sin(2πkn/N)取值0、±1的點(diǎn)，如圖中x點(diǎn)復(fù)共軛對(duì)稱性(圖中o點(diǎn))： WNk(N-n)=WN-kn=(WNkn)*對(duì)n,k的周期性： WNkn=WNk(N+n)=WN(N+k)n運(yùn)算量∝N2，當(dāng)N=PQ時(shí)，可分解為P個(gè)Q點(diǎn)DFT59.2Goertzel算法∵WN-kN=ej(2π/N)kN=1，且當(dāng)n在[0,N)外時(shí)x[n]=0∴X[k]=x[r]WNkr(WN-kNu[N-r])=x[r]WN-k(N-r)u[N-r]=yk[n]|n=N式中：yk[n]=x[r]WN-k(n-r)u[n-r]=x[n]*(WN-knu[n])Hk(z)=1/(1-WN-kz-1)=(1-WNkz-1)/(1-2cos(2πk/N)z-1+z-2)yk[n]是x[n]經(jīng)Hk(z)的響應(yīng)，X[k]是其n=N時(shí)的輸出∑r=0N-1∑r=-∞∞∑r=-∞∞6Goertzel算法信號(hào)流圖實(shí)現(xiàn)初始松弛條件等效有限長(zhǎng)x[n]=0，n<0復(fù)數(shù)x[n]時(shí)每次遞推要4次實(shí)乘、4次實(shí)加為求X[k]=yk[N]，需遞推N次，4N實(shí)乘、4N實(shí)加計(jì)算全部X[k](k=0,…,N-1)時(shí)運(yùn)算量并未減少優(yōu)點(diǎn)：不需要計(jì)算或存儲(chǔ)WNkn易于計(jì)算單個(gè)或幾個(gè)X[k]單反饋回路，結(jié)構(gòu)簡(jiǎn)單7Goertzel算法信號(hào)流圖改進(jìn)反饋迭代N次后，僅最后需算一次零點(diǎn)回路極點(diǎn)回路僅一個(gè)實(shí)系數(shù)，實(shí)乘次數(shù)減少近半X[N-k]運(yùn)算類同X[k]，僅零點(diǎn)系數(shù)為(-WNk)*單X[k]運(yùn)算量：2N+4次實(shí)乘4N+4次實(shí)加N個(gè)X[k]運(yùn)算量：N(N+4)次實(shí)乘N(2N+4)次實(shí)加8Goertzel算法特點(diǎn)系數(shù)存儲(chǔ)僅一復(fù)一實(shí)大為減少，迭代中隱含算出結(jié)構(gòu)簡(jiǎn)單，資源開銷小，非常適合硬件實(shí)現(xiàn)與DFT直接計(jì)算相比，單點(diǎn)X[k]運(yùn)算量減少近兩倍，總運(yùn)算量減少近4倍與基2FFT相比，N為任意值不受限由于X[k]可單點(diǎn)計(jì)算，當(dāng)僅需求出M個(gè)X[k]，且M<log2N時(shí)，運(yùn)算量小于基2FFT99.3按時(shí)間抽取的FFT算法基2為例：N=2v，可將DFT輸入按偶奇對(duì)分X[k]=x[n]WNkn=x[2r]WN2rk+x[2r+1]WN(2r+1)k=x[2r]WN/2rk+WNkx[2r+1]WN/2rk=G[k]+WNkH[k]注意：∵WNm=e-jm2π/N=WN/m，∴WN2=WN/2

N點(diǎn)DFTX[k]可分解為兩個(gè)N/2點(diǎn)的DFTG[k]與H[k]的線性組合G[k]、H[k]的周期為N/2，G[k]=G[k-N/2]∑r=0N/2-1∑r=0N/2-1∑r=0N/2-1∑r=0N/2-110N點(diǎn)DFT分解為兩個(gè)N/2點(diǎn)DFT11N/2點(diǎn)DFT可分解為N/4點(diǎn)DFT由于基2，N/2點(diǎn)DFTG[k]、H[k]可繼續(xù)按輸入的“偶奇”分解為各兩個(gè)N/4點(diǎn)DFT注意輸入的序號(hào)、系數(shù)的下標(biāo)與冪次12兩級(jí)分解后的完整流圖13繼續(xù)對(duì)分，直至最小兩點(diǎn)DFT基2時(shí)，有v=log2N級(jí)系數(shù)標(biāo)注規(guī)律2點(diǎn)DFT：W20,W214點(diǎn)DFT：W40,W41,W42,W43N點(diǎn)DFT：WN0,WN1,…,WNN-1為使用統(tǒng)一系數(shù)表，將下標(biāo)都調(diào)整為N：Wmk=WNkN/m14一個(gè)8點(diǎn)DFT完全分解后流圖15基本蝶形單元的優(yōu)化整個(gè)流圖由左下圖基本蝶形單元構(gòu)成，故基2算法又常稱為蝶形算法，蝶形流圖每個(gè)蝶形單元的兩個(gè)系數(shù)的冪一定相差N/2：∵WNN/2=e-j2π/2=-1，∴WNr+N/2=WNrWNN/2=-WNr乘法可合并，優(yōu)化蝶形單元如右下圖16輸出正位序同址計(jì)算FFT流圖17同址計(jì)算按優(yōu)化蝶形單元，兩個(gè)輸入、兩個(gè)輸出作為一組計(jì)算，結(jié)果存放回原址流圖表現(xiàn)為水平蝶形優(yōu)點(diǎn)是對(duì)存儲(chǔ)要求最小，N點(diǎn)FFT僅需一組N個(gè)復(fù)變量存儲(chǔ)單元(可能另需N/2個(gè)復(fù)系數(shù)存儲(chǔ)單元，復(fù)數(shù)一般按實(shí)部、虛部分別存放，非同址一般需二組變量存儲(chǔ)單元)18基2FFT的運(yùn)算量基2時(shí)，分解級(jí)數(shù)為v=log2N級(jí)，且每級(jí)N/2次復(fù)乘，N次復(fù)加N點(diǎn)DFT總運(yùn)算量：(N/2)log2N復(fù)乘，Nlog2N復(fù)加與DFT直接計(jì)算∝N2相比，N較大時(shí)效果極佳當(dāng)利用WNN=WN0=1后，復(fù)乘還可減少。但采用同址計(jì)算編程實(shí)現(xiàn)時(shí)，為保持過程統(tǒng)一，減少程序量，有時(shí)未進(jìn)行乘加的極致減少19倒位序?qū)ぶ放c偶奇對(duì)分過程將序列的序號(hào)n用2進(jìn)制表達(dá)為n2n1n0第一次進(jìn)行最后一級(jí)偶奇分解實(shí)際就是按最低位n0排序第二次按n1排序…規(guī)律:輸入x[n2n1n0]存儲(chǔ)在第[n0n1n2]行20FFT流圖的拓?fù)渥冃位?FFT流圖屬多輸入、多輸出信號(hào)流圖只要流圖的節(jié)點(diǎn)關(guān)系、傳輸比不變，代表的運(yùn)算關(guān)系就不會(huì)改變拓?fù)鋵W(xué)的變形有許多種，影響的是數(shù)據(jù)的讀取與存儲(chǔ)次序，以及硬件實(shí)現(xiàn)或軟件編程的架構(gòu)最常用的還是輸入或輸出正位序的同址計(jì)算FFT信號(hào)流圖算法21按時(shí)間抽取輸入正位序FFT22例：輸入輸出皆正位序的流圖中間過程復(fù)雜，無(wú)實(shí)用意義23例：每級(jí)同架構(gòu)的FFT流圖支持?jǐn)?shù)據(jù)順序讀寫，可用于超長(zhǎng)序列DFT計(jì)算249.4按頻率抽取的FFT算法同樣基2為例，可將DFT輸出取偶數(shù)部分X[2r]=x[n]WN2rn=x[n]WN2rn+x[n]WN2rn={x[n]WN2rn+x[n+N/2]WN2r(n+N/2)}={x[n]+x[n+N/2]}WN/2rn類似得奇數(shù)部分：X[2r+1]={x[n]-x[n-N/2]}WNnWN/2rn輸出減抽樣，導(dǎo)致輸入信號(hào)混疊∑n=0N/2-1∑n=N/2N-1∑n=0N/2-1∑n=0N/2-1∑n=0N/2-125按輸出將N點(diǎn)DFT分解為二N/2點(diǎn)26輸出逐次二分，直至最小二點(diǎn)27按頻率抽取輸入正位序FFT28按頻率抽取FFT特點(diǎn)運(yùn)算量與按時(shí)間抽取FFT完全相同：(N/2)log2N復(fù)乘，Nlog2N復(fù)加按時(shí)間與按頻率抽取FFT間關(guān)系：多輸入/輸出信號(hào)流圖的整體倒置形式上的流圖倒置，非單輸入/單輸出的信號(hào)流圖倒置定理每種按時(shí)間抽取的FFT流圖，必有對(duì)應(yīng)倒置的按頻率抽取的FFT流圖節(jié)點(diǎn)關(guān)系、系數(shù)標(biāo)注完全對(duì)稱，規(guī)律相似29例：流圖的整體倒置按時(shí)間抽取輸出正位序按頻率抽取輸入正位序309.5FFT實(shí)現(xiàn)問題討論IDFT的差別：系數(shù)WN-i=(WNi)*，復(fù)數(shù)運(yùn)算虛部反號(hào)，可與DFT用同一系數(shù)表；結(jié)果差1/N因子，僅僅為定標(biāo)問題倒位序的實(shí)現(xiàn)：DSP處理器中有專門的硬件，效率極高系數(shù)可遞推算出，更常用表格，包括倒位序也可預(yù)安排N非基2時(shí)，常用補(bǔ)零至基2形式復(fù)合數(shù)N=PQ也可導(dǎo)出類似算法31同址計(jì)算FFT流圖規(guī)律輸入正位序，輸出一定倒位序；反之亦然級(jí)數(shù)m=1,2,…,log2N輸入倒位序時(shí)，第m級(jí)數(shù)據(jù)間隔為2m-1(即：1,2,4,8,…,N/2)輸入正位序時(shí)，第m級(jí)數(shù)據(jù)間隔為N/2m按時(shí)間抽取，每級(jí)系數(shù)為W2m=WNN/2m的冪次，位序同輸出按頻率抽取，注意倒置關(guān)系，系數(shù)冪次位序同輸入32例：用FFT實(shí)現(xiàn)塊卷積DFT采用按頻率抽取輸入正位序IDFT采用按時(shí)間抽取輸出正位序二者系數(shù)冪次也都正位序，統(tǒng)一且方便運(yùn)算總的輸入/輸出都正位序不必調(diào)整h[n]預(yù)先用DFT求出H[k]保存頻域雖然是倒位序但位序相同，全部對(duì)應(yīng)點(diǎn)相乘可不管位序339.6有限寄存器長(zhǎng)度的影響采用簡(jiǎn)化線性模型分析有限字長(zhǎng)定點(diǎn)運(yùn)算誤差產(chǎn)生的影響及優(yōu)化方法分析以基2按時(shí)間抽取FFT算法為例將運(yùn)算誤差作為加性噪聲處理如下圖小數(shù)格式定點(diǎn)運(yùn)算字長(zhǎng)：B+1位34單復(fù)乘節(jié)點(diǎn)的噪聲方差乘法舍入近似為統(tǒng)計(jì)獨(dú)立的誤差噪聲設(shè)噪聲在區(qū)間(-2-B/2,2-B/2)均勻分布復(fù)乘對(duì)應(yīng)四個(gè)實(shí)乘，誤差噪聲相互無(wú)關(guān)所有的誤差噪聲與輸入/輸出信號(hào)不相關(guān)單復(fù)乘節(jié)點(diǎn)噪聲方差：E{|e|2}=4x2-2B/12=2-2B/3=σB235任一輸出的噪聲疊加路徑36任一單輸出值的總噪聲任一輸出都有(N-1)個(gè)相關(guān)的復(fù)乘節(jié)點(diǎn)每個(gè)復(fù)乘節(jié)點(diǎn)噪聲特性相同，統(tǒng)計(jì)無(wú)關(guān)(未考慮WN0=1的節(jié)點(diǎn)其實(shí)并無(wú)誤差)|WNr|=1，前級(jí)噪聲通過后級(jí)回路不改變方差任一單輸出值總噪聲方差：E{|F(k)|2}=(N-1)σB2≈

NσB237運(yùn)算無(wú)溢出分析注意|WNr|=1，由FFT基本蝶形單元可得：max(|Xm-1[p]|,|Xm-1[q]|)≤max(|Xm[p]|,|Xm[q]|)2max(|Xm-1[p]|,|Xm-1[q]|)≥max(|Xm[p]|,|Xm[q]|)幅度逐級(jí)遞增，輸出無(wú)溢出，中間定無(wú)溢出運(yùn)算無(wú)溢出|X[k]|<1的充分條件：|x[n]|<1/N，0≤n,k≤N-1Xm[p]Xm[q]Xm-1[p]Xm-1[q]-Xm-1[q]38無(wú)溢出時(shí)單輸出信號(hào)均方幅度設(shè)輸入信號(hào)x[n]為白噪聲復(fù)序列滿足運(yùn)算無(wú)溢出：Re[x[n]]與Im[x[n]]<2-1/2/N設(shè)虛、實(shí)部在區(qū)間(-2-1/2/N,2-1/2/N)均勻分布輸入x[n]均方幅度：E{|x[n]|2}=1/3N2=σx2運(yùn)算無(wú)溢出時(shí)單輸出信號(hào)X[k]均方幅度：E{|X[k]|2}=Nσx2=1/3N39FFT定點(diǎn)運(yùn)算信噪比分析信噪比：SNR=E{|X[k]|2}/E{|F(k)|2}SNR=(1/3N)/NσB2=22B/N2=22(B-V)SNR反比于N2，N擴(kuò)大一倍，SNR將下降四倍SNR正比于4B，B增加一