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

CH9DFT計算實現(xiàn)方法—FFT9.1DFT運算量估計及算法分析9.2Goertzel算法9.3按時間抽取的FFT算法9.4按頻率抽取的FFT算法9.5FFT實現(xiàn)問題討論9.6有限寄存器長度的影響19.1DFT運算量估計及算法分析算法優(yōu)劣的衡量與應用相關，工程中經(jīng)?？紤]的有：精度、信噪比實時性、收斂性存儲量、硬件成本運算量：乘法與加法次數(shù)各項要求有時矛盾，研究人員在實現(xiàn)中要做好和諧工作2例：三角函數(shù)的實現(xiàn)實時性要求很高時用查表法，如雷達PPI極坐標到光柵掃描直角坐標的轉換通用性要求較好時常用冪級數(shù)展開法，如一般計算機程序的函數(shù)庫輸出有規(guī)則時可用遞推法，如波形發(fā)生器運算量實時性存儲量靈活性查表法小極好很大一般Z變換遞推法中很好小差冪級數(shù)展開法大較好較小好3DFT運算量估計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復數(shù)N點DFT運算：N2次復乘，N(N-1)次復加實數(shù)N點DFT運算：4N2次實乘，N(4N-2)次實加4DFT算法分析算法分析的目的是減少乘法與加法的次數(shù)Re[WNkn]=cos(2πkn/N)，Im[WNkn]=-sin(2πkn/N)取值0、±1的點，如圖中x點復共軛對稱性(圖中o點)： WNk(N-n)=WN-kn=(WNkn)*對n,k的周期性： WNkn=WNk(N+n)=WN(N+k)n運算量∝N2，當N=PQ時，可分解為P個Q點DFT59.2Goertzel算法∵WN-kN=ej(2π/N)kN=1，且當n在[0,N)外時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)的響應，X[k]是其n=N時的輸出∑r=0N-1∑r=-∞∞∑r=-∞∞6Goertzel算法信號流圖實現(xiàn)初始松弛條件等效有限長x[n]=0，n<0復數(shù)x[n]時每次遞推要4次實乘、4次實加為求X[k]=yk[N]，需遞推N次，4N實乘、4N實加計算全部X[k](k=0,…,N-1)時運算量并未減少優(yōu)點：不需要計算或存儲WNkn易于計算單個或幾個X[k]單反饋回路，結構簡單7Goertzel算法信號流圖改進反饋迭代N次后，僅最后需算一次零點回路極點回路僅一個實系數(shù)，實乘次數(shù)減少近半X[N-k]運算類同X[k]，僅零點系數(shù)為(-WNk)*單X[k]運算量：2N+4次實乘4N+4次實加N個X[k]運算量：N(N+4)次實乘N(2N+4)次實加8Goertzel算法特點系數(shù)存儲僅一復一實大為減少，迭代中隱含算出結構簡單，資源開銷小，非常適合硬件實現(xiàn)與DFT直接計算相比，單點X[k]運算量減少近兩倍，總運算量減少近4倍與基2FFT相比，N為任意值不受限由于X[k]可單點計算，當僅需求出M個X[k]，且M<log2N時，運算量小于基2FFT99.3按時間抽取的FFT算法基2為例：N=2v，可將DFT輸入按偶奇對分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點DFTX[k]可分解為兩個N/2點的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點DFT分解為兩個N/2點DFT11N/2點DFT可分解為N/4點DFT由于基2，N/2點DFTG[k]、H[k]可繼續(xù)按輸入的“偶奇”分解為各兩個N/4點DFT注意輸入的序號、系數(shù)的下標與冪次12兩級分解后的完整流圖13繼續(xù)對分，直至最小兩點DFT基2時，有v=log2N級系數(shù)標注規(guī)律2點DFT：W20,W214點DFT：W40,W41,W42,W43N點DFT：WN0,WN1,…,WNN-1為使用統(tǒng)一系數(shù)表，將下標都調(diào)整為N：Wmk=WNkN/m14一個8點DFT完全分解后流圖15基本蝶形單元的優(yōu)化整個流圖由左下圖基本蝶形單元構成，故基2算法又常稱為蝶形算法，蝶形流圖每個蝶形單元的兩個系數(shù)的冪一定相差N/2：∵WNN/2=e-j2π/2=-1，∴WNr+N/2=WNrWNN/2=-WNr乘法可合并，優(yōu)化蝶形單元如右下圖16輸出正位序同址計算FFT流圖17同址計算按優(yōu)化蝶形單元，兩個輸入、兩個輸出作為一組計算，結果存放回原址流圖表現(xiàn)為水平蝶形優(yōu)點是對存儲要求最小，N點FFT僅需一組N個復變量存儲單元(可能另需N/2個復系數(shù)存儲單元，復數(shù)一般按實部、虛部分別存放，非同址一般需二組變量存儲單元)18基2FFT的運算量基2時，分解級數(shù)為v=log2N級，且每級N/2次復乘，N次復加N點DFT總運算量：(N/2)log2N復乘，Nlog2N復加與DFT直接計算∝N2相比，N較大時效果極佳當利用WNN=WN0=1后，復乘還可減少。但采用同址計算編程實現(xiàn)時，為保持過程統(tǒng)一，減少程序量，有時未進行乘加的極致減少19倒位序?qū)ぶ放c偶奇對分過程將序列的序號n用2進制表達為n2n1n0第一次進行最后一級偶奇分解實際就是按最低位n0排序第二次按n1排序…規(guī)律:輸入x[n2n1n0]存儲在第[n0n1n2]行20FFT流圖的拓撲變形基2FFT流圖屬多輸入、多輸出信號流圖只要流圖的節(jié)點關系、傳輸比不變，代表的運算關系就不會改變拓撲學的變形有許多種，影響的是數(shù)據(jù)的讀取與存儲次序，以及硬件實現(xiàn)或軟件編程的架構最常用的還是輸入或輸出正位序的同址計算FFT信號流圖算法21按時間抽取輸入正位序FFT22例：輸入輸出皆正位序的流圖中間過程復雜，無實用意義23例：每級同架構的FFT流圖支持數(shù)據(jù)順序讀寫，可用于超長序列DFT計算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輸出減抽樣，導致輸入信號混疊∑n=0N/2-1∑n=N/2N-1∑n=0N/2-1∑n=0N/2-1∑n=0N/2-125按輸出將N點DFT分解為二N/2點26輸出逐次二分，直至最小二點27按頻率抽取輸入正位序FFT28按頻率抽取FFT特點運算量與按時間抽取FFT完全相同：(N/2)log2N復乘，Nlog2N復加按時間與按頻率抽取FFT間關系：多輸入/輸出信號流圖的整體倒置形式上的流圖倒置，非單輸入/單輸出的信號流圖倒置定理每種按時間抽取的FFT流圖，必有對應倒置的按頻率抽取的FFT流圖節(jié)點關系、系數(shù)標注完全對稱，規(guī)律相似29例：流圖的整體倒置按時間抽取輸出正位序按頻率抽取輸入正位序309.5FFT實現(xiàn)問題討論IDFT的差別：系數(shù)WN-i=(WNi)*，復數(shù)運算虛部反號，可與DFT用同一系數(shù)表；結果差1/N因子，僅僅為定標問題倒位序的實現(xiàn)：DSP處理器中有專門的硬件，效率極高系數(shù)可遞推算出，更常用表格，包括倒位序也可預安排N非基2時，常用補零至基2形式復合數(shù)N=PQ也可導出類似算法31同址計算FFT流圖規(guī)律輸入正位序，輸出一定倒位序；反之亦然級數(shù)m=1,2,…,log2N輸入倒位序時，第m級數(shù)據(jù)間隔為2m-1(即：1,2,4,8,…,N/2)輸入正位序時，第m級數(shù)據(jù)間隔為N/2m按時間抽取，每級系數(shù)為W2m=WNN/2m的冪次，位序同輸出按頻率抽取，注意倒置關系，系數(shù)冪次位序同輸入32例：用FFT實現(xiàn)塊卷積DFT采用按頻率抽取輸入正位序IDFT采用按時間抽取輸出正位序二者系數(shù)冪次也都正位序，統(tǒng)一且方便運算總的輸入/輸出都正位序不必調(diào)整h[n]預先用DFT求出H[k]保存頻域雖然是倒位序但位序相同，全部對應點相乘可不管位序339.6有限寄存器長度的影響采用簡化線性模型分析有限字長定點運算誤差產(chǎn)生的影響及優(yōu)化方法分析以基2按時間抽取FFT算法為例將運算誤差作為加性噪聲處理如下圖小數(shù)格式定點運算字長：B+1位34單復乘節(jié)點的噪聲方差乘法舍入近似為統(tǒng)計獨立的誤差噪聲設噪聲在區(qū)間(-2-B/2,2-B/2)均勻分布復乘對應四個實乘，誤差噪聲相互無關所有的誤差噪聲與輸入/輸出信號不相關單復乘節(jié)點噪聲方差：E{|e|2}=4x2-2B/12=2-2B/3=σB235任一輸出的噪聲疊加路徑36任一單輸出值的總噪聲任一輸出都有(N-1)個相關的復乘節(jié)點每個復乘節(jié)點噪聲特性相同，統(tǒng)計無關(未考慮WN0=1的節(jié)點其實并無誤差)|WNr|=1，前級噪聲通過后級回路不改變方差任一單輸出值總噪聲方差：E{|F(k)|2}=(N-1)σB2≈

NσB237運算無溢出分析注意|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]|)幅度逐級遞增，輸出無溢出，中間定無溢出運算無溢出|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無溢出時單輸出信號均方幅度設輸入信號x[n]為白噪聲復序列滿足運算無溢出：Re[x[n]]與Im[x[n]]<2-1/2/N設虛、實部在區(qū)間(-2-1/2/N,2-1/2/N)均勻分布輸入x[n]均方幅度：E{|x[n]|2}=1/3N2=σx2運算無溢出時單輸出信號X[k]均方幅度：E{|X[k]|2}=Nσx2=1/3N39FFT定點運算信噪比分析信噪比：SNR=E{|X[k]|2}/E{|F(k)|2}SNR=(1/3N)/NσB2=22B/N2=22(B-V)SNR反比于N2，N擴大一倍，SNR將下降四倍SNR正比于4B，B增加一