綜合課程設(shè)計報告第三組

綜合課程設(shè)計報告

基于ADSP-2191的數(shù)字信號處理實現(xiàn)題目：基二頻率抽取的快速傅立葉變換算法小組成員:安海鵬04008220高軒04008213 王康蔚04008215朱峻崗04008212蔣駿寧（組長）04008236東南大學信息科學與工程學院二○一一年一月綜合課程設(shè)計—基于ADSP-2191的數(shù)字信號處理實現(xiàn)-PAGE4-題目概述快速傅立葉變換是一種快速有效地計算離散傅立葉變換的方法。一般情況下輸入的時間序列及其離散傅立葉變換是用復數(shù)表示的，直接計算DFT及IDFT需要次復數(shù)乘法以及次復數(shù)加法，由于1次復數(shù)乘法對應4次實數(shù)懲罰和2次實數(shù)加法，1次復數(shù)加法需要做2次實數(shù)加法，所以1次離散傅立葉變換需要做次實數(shù)乘法及次實數(shù)加法。20世紀60年代Cooley和Tukey提出了一種離散傅立葉變換的快速算法，將其運算量下降為次復數(shù)乘法和次復數(shù)加法。本次基于ADSP-2191的數(shù)字信號處理實現(xiàn)基二頻率抽取的算法就是對上述功能在DSP平臺上的利用匯編語言的實現(xiàn)。算法介紹基二頻率抽取FFT運算流圖圖2.1.1圖2.1.1為N=16的基二頻率抽取的FFT運算流圖，構(gòu)成了課題目標算法的理論基礎(chǔ)，流圖包含：蝶形運算和序數(shù)重排兩個部分，總的蝶形個數(shù)為，果重新排序得到正確的結(jié)果。在N=16的情況下，一共有32次蝶形運算，即一共需要32次復數(shù)乘法運算和64次復數(shù)加法運算。圖2.1.2輸入稀疏度為4的N=16的基二頻率抽取的FFT運算流圖圖2.1.2為本課程要求的輸入稀疏度為4基二頻率抽取的FFT運算流圖，從圖上可以得到需要進行12次半蝶形運算和16次蝶形運算，一共需要進行28次復數(shù)乘法和32次復數(shù)加法運算。對比圖2.1.1的流圖，N=16的FFT運算流圖需要32次復數(shù)乘法和64次復數(shù)加法，圖2.蝶形運算與半蝶形運算復數(shù)加法與減法復數(shù)加法一次復數(shù)加法包含行兩次實數(shù)加法過程。對于實數(shù)加法，設(shè)輸入X操作數(shù)X,Y操作數(shù)Y，X,Y，且X,Y均為16位有符號數(shù)，加法得到的結(jié)果Z,Z為17位有符號數(shù)。如果輸入的X,Y操作數(shù)均為，得到的結(jié)果就是溢出的，如果將存入寄存器繼續(xù)運算的話，將帶來鏈式錯誤?？紤]上述情況，需要Z和X,Y具有相同的范圍即，本程序采取先將X,Y均算術(shù)右移一位，即可使得Z的范圍滿足要求。同時考慮到，最后一位的舍入誤差，在右移后，先將X,Y的最后一位進行加法運算，在考慮進位的情況下，進行主值加法。此方法的誤差為不考慮末位進位的1/2，降低的誤差率。復數(shù)減法一次復數(shù)減法包含兩次實數(shù)減法過程。對于實數(shù)減法，設(shè)輸入X操作數(shù)X,Y操作數(shù)Y，X,Y，且X,Y均為16位有符號數(shù)，加法得到的結(jié)果Z,Z為17位有符號數(shù)。如果輸入的X,Y操作數(shù)為和，得到的結(jié)果就是溢出的，如果將存入寄存器繼續(xù)運算的話，將帶來鏈式錯誤?？紤]上述情況，需要Z和X,Y具有相同的范圍即，本程序采取先將X,Y均算術(shù)右移一位，即可使得Z的范圍滿足要求。同時考慮到，最后一位的舍入誤差，在右移后，先將X,Y的最后一位進行減法運算，在考慮借位的情況下，進行主值減法。此方法的誤差為不考慮末位借位的1/2，降低的誤差率。復數(shù)乘法一次復數(shù)乘法包含四次實數(shù)乘法和二次實數(shù)加法，設(shè)輸入數(shù)據(jù)，為一次復數(shù)乘法運算的結(jié)果，其中為的實部和虛部，其絕對值都不超過1。為輸入數(shù)據(jù)的實部和虛部，其均的范圍。在極端情況下，假設(shè)均為，在阻止溢出的情況下，可以采取復數(shù)加減法的策略，即對進行算術(shù)右移后再進行運算，也可以預先設(shè)計的輸入數(shù)格式，將通用的1.15的表示格式變成2.14的表示格式。這樣，在輸入均為的情況下，依然能夠防止輸出溢出。序數(shù)重排表2.3.1N=16基二頻率抽取序號重排位置序號與實際序號輸出數(shù)據(jù)位置序號輸出數(shù)據(jù)實際序號0000X(0)00000001X(8)10000010X(4)01000011X(12)11000100X(2)00100101X(10)10100110X(6)01100111X(14)11101000X(1)00011001X(9)10011010X(5)01011011X(13)11011100X(3)00111101X(11)10111110X(7)01111111X(15)1111模塊劃分整個程序分為半蝶形運算模塊、蝶形運算模塊、比特反轉(zhuǎn)模塊、源地址與目標地址交換模塊。各個模塊包含的運算如下表：表3.1模塊功能表實數(shù)加法實數(shù)減法復數(shù)乘法數(shù)據(jù)讀寫寄存器讀寫半蝶形運算模塊無無有有有蝶形運算模塊有有有有有比特反轉(zhuǎn)模塊無無無有有源地址與目標地址交換模塊無無無無有地址重置與操作間隔改寫模塊無無無無有實數(shù)加法：為2.2.1.1復數(shù)加法中描述的實數(shù)加法。實數(shù)減法：為2.2.1.2復數(shù)減法中描述的實數(shù)減法。復數(shù)乘法：為2.2.2復數(shù)乘法中描述的復數(shù)乘法，因ADSP2191的乘法加法累加器（MAC）特性，故實數(shù)乘法中包含的復數(shù)加法減法與2.2.1.1與2.2.1.2中描述無關(guān)。數(shù)據(jù)讀寫：利用ADSP2191尋址方式，涉及數(shù)據(jù)寄存器和存儲器之間的數(shù)據(jù)交互。寄存器操讀寫：包含寄存器讀入寫出，寄存器特殊位賦值操作。半蝶形運算模塊圖3.1.1一組N=4的半蝶形運算一次半蝶形運算包含一次數(shù)據(jù)移動和一次復數(shù)乘法運算。全蝶形運算模塊圖3.2.1一組N=2的全蝶形運算一次全蝶形運算包含一次復數(shù)加法、一次復數(shù)減法、一次復數(shù)乘法。比特反轉(zhuǎn)模塊利用ADSP2191自帶的比特反轉(zhuǎn)功能實現(xiàn)2.3中的序數(shù)重排功能。源地址與目標地址交換模塊為了提高ADSP系統(tǒng)存儲器的利用率，將上次半蝶形或蝶形運算的源地址作為本次半蝶形或蝶形運算的目標地址，可以減少存儲單元的浪費。地址重置與操作間隔改寫模塊由于中的數(shù)據(jù)在每次運算中由固定地址讀入，而非將數(shù)據(jù)全部存儲于寄存器中，每次蝶形操作或半蝶形操作后，見圖2.1.1需要將數(shù)據(jù)指針重新置于頂部，叫做地址重置。操作間隔改寫在于半蝶型運算中一個輸入數(shù)據(jù)需要同時作為兩個目標數(shù)據(jù)的輸入數(shù)據(jù)，而且兩個目標數(shù)據(jù)的間隔在2級半蝶形運算間是變化的。同時，2級全蝶形運算中，一次加法和一次減法乘法的輸入數(shù)據(jù)間隔和目標數(shù)據(jù)間隔也不一樣，所以需要改寫操作間隔。定點數(shù)處理和各模塊中定點數(shù)的配置對于，如2.2.1.2所述，將其定義為2.14格式。對于4個輸入數(shù)據(jù)，其實部和虛部均為4.12。一次復數(shù)加法、減法、乘法都將提高數(shù)據(jù)的位數(shù)，即4.12格式成為5.11格式，5.11格式成為6.10格式，以此類推。圖4.1塊浮點數(shù)下的數(shù)據(jù)格式上圖為在極端情況下，每次復數(shù)加法、減法、乘法都存在進位的情況下的塊浮點數(shù)格式，輸出的數(shù)據(jù)從6.10格式到10.06格式，在理論精度達到最大的情況下，帶來了格式的不統(tǒng)一，數(shù)據(jù)操作的困難。圖4.2定點數(shù)下數(shù)據(jù)格式上圖為在極端情況下，每次復數(shù)加法、減法、乘法都存在進位的情況下的定點數(shù)格式，每一級輸出的格式都是相同的，最后的輸出格式為10.06。定點實現(xiàn)與浮點結(jié)果對比測試序列：（4.12輸入格式）I路（實數(shù)部分）：Q路（虛數(shù)部分）：用復數(shù)表示為：塊浮點計算結(jié)果與MATLAB仿真對比仿真I計算I(MATLAB)誤差仿真Q計算Q(MATLAB)誤差FC00FC000FC00FC0000CF70CF8029302941B361B37F4BAF4BA00B0D0B0D0EF74EF75000000000D800D8000F7EAF7EA0F3E6F3E60FB77FB770FE1FFE20030A030BFF65FF66160016000EA00EA000040004000EBE1EBE10F865F8650EDA3EDA30FA0FFA0F0FC07FC08FA00FA00000000000000CD00CD0FF77FF7801440145FCF0FCF1FFAFFFB0FE6CFE6DMALTAB計算使用ROUND函數(shù)，CEIL,FLOOR差距均不大參考流圖4.1中所對應的數(shù)據(jù)格式，所有的誤差均不超過對應塊浮點最后一位的絕對值，誤差來源有：精度誤差，即值用2.14格式的16位數(shù)表示帶來的誤差；舍入誤差，因為運算進位的關(guān)系，若干位值被自動忽略。定點計算結(jié)果與MATLAB仿真對比仿真I計算I(MATLAB)誤差仿真Q計算Q(MATLAB)誤差FFC0FFC00FFC0FFC00019E019F0052005303660367FE97FF97002C302C4FBDDFBDD0000000000FB00FB000FDFAFDFA0FCF9FCF90FEDDFEDEFF87FF88018501850FFB2FFB302C002C00FD40FD400010001000FAF8FAF80FE19FE190FB68FB69FD07FD070FE04FE040FE80FE800000000000006600660FFBBFFBC00A200A20FE78FE780FFAFFFB0FE6CFE6C0MALTAB計算使用ROUND函數(shù)，CEIL,FLOOR差距均不大。參考流圖4.2中所對應的數(shù)據(jù)格式，所有的誤差均不超過對應塊浮點最后一位的絕對值，誤差來源有：精度誤差，即值用2.14格式的16位數(shù)表示帶來的誤差；舍入誤差，因為運算進位的關(guān)系，若干位值被自動忽略。算法復雜度統(tǒng)計對于最終的程序，即定點數(shù)法圖4.2流圖，經(jīng)過統(tǒng)計一共進行了16次復數(shù)加法，16次復數(shù)減法，28次復數(shù)乘法，28次位移操作。一共進行了3次半蝶形函數(shù)調(diào)用，每次進行4次半蝶形運算；8次全蝶形函數(shù)調(diào)用每次進行2次全蝶形運算；1次比特反轉(zhuǎn)。1次實數(shù)加法對應8個指令周期，1次實數(shù)減法對應8個指令周期，1次復數(shù)乘法對應5個指令周期。主體含103個指令周期（不含最后RTS），比特反轉(zhuǎn)占74個指令周期，一次半蝶形函數(shù)含91個指令周期，一個全蝶形函數(shù)含135個指令周期，整個程序調(diào)用需要1530個指令周期。小結(jié)這次基于ADSP2191平臺的課程設(shè)計給予了我們一次用實踐檢驗理論正確性的機會，增進了團隊間的合作，提高了小組成員間協(xié)作分工的能力。在此過程中，實際中的很多限制，比如定點數(shù)等給予了我們在設(shè)計算法中提出了許多挑戰(zhàn)，也讓增進了我們，在現(xiàn)實環(huán)境中去應用理論的能力參考文獻[1]吳鎮(zhèn)揚數(shù)字信號處理第二版高等教育出版社[2]胡廣書數(shù)字信號處理理論算法與實現(xiàn)第二版清華大學出版社源代碼/************************************************************************//filename:test.asm//purpose:bodyofprogram#definenum4#definecoe8#definegap1=8;#definegap2=4;#definegap3=2;.externadd;.externsub;.externmulti;.externhalfbutterfly;.externtotalbutterfly;.externsoudes;.externbitrev;/*DMdata*/.section/datadata1;.varinput_i[num]="input_i.dat"; //inputitakes8000-8003.varinput_q[num]="input_q.dat";//inputqtakes8004-8007.var coeff_i[coe]="coeff_i.dat";//coeffitakes8008-800f.var coeff_q[coe]="coeff_q.dat";//coeffqtakes8010-8017/*PMinterruptvectorcode*/.section/pmIVreset;JUMPstart;NOP;NOP;NOP;/*Programmemorycode*/.section/pmprogram;start: i0=input_i; i1=input_q; i2=coeff_i; i3=coeff_q; m0=1; // m3=1;//forcos m4=1;//increaseindatareadandwrite m5=8;//gap m6=1;//increasebetweentwo i4=0x8040;//resultaddressfori i5=0x8050; //resultaddressforq //layerone callhalfbutterfly; //atthistimewedon'tneedi0i1sowechoosei0i1asournextdestinationaddress //i0=0x8040; //i1=0x8050; //i4=0x8020; //i5=0x8030; //callsoudes; //callcossinrst; i2=coeff_i;//returncostostart i3=coeff_q;//returnsintostart //layertwo i0=0x8040; i1=0x8050; i4=0x8020; i5=0x8030; m3=2; m5=4;callhalfbutterfly;// callcossinrst; i2=coeff_i;//returncostostart i3=coeff_q;//returnsintostart ax0=dm(i0+=4); ax0=dm(i1+=4); ax0=dm(i4+=4); ax0=dm(i5+=4); callhalfbutterfly; //callcossinrst; i2=coeff_i;//returncostostart i3=coeff_q;//returnsintostart //callsoudes; i0=0x8020; i1=0x8030; i4=0x8040; i5=0x8050; //layerthree m0=1;//gapbetweentwooutput m1=2;//gapininput m5=2;//gapinoutput m6=1;//gapbetweentwooutput m3=4; cntr=4; docir_3untilce; calltotalbutterfly; ax0=dm(i0+=2); ax0=dm(i1+=2); ax0=dm(i4+=2); ax0=dm(i5+=2); // callcossinrst; i2=coeff_i;//returncostostart i3=coeff_q;//returnsintostart cir_3:nop; // callsoudes; i0=0x8040; i1=0x8050; i4=0x8020; i5=0x8030; //layerfour m3=0; m0=2; m1=1; m5=1; m6=2; cntr=4; docir_4untilce; calltotalbutterfly; // callcossinrst; i2=coeff_i;//returncostostart i3=coeff_q;//returnsintostart //ax0=dm(i0+=2) cir_4:nop; //callsoudes; i4=0x8020; i5=0x8030; i0=0x0201; i1=0x0a01; //bitrev callbitrev;rts; //filename:add.asm//purpose:addopcode//.globaladd;.globaladd;/*programmemorycode*///input:ax0,ay0,outputar.section/pmprogram;add: se=-1; sr=ashiftax0(hi); ax0=sr1;ax1=sr0; sr=ashiftay0(hi); ay0=sr1;ay1=sr0; ar=ax1+ay1; ar=ax0+ay0+c; rts; //filename:sub.asm//purpose:subopcode//inputax0,ay0,outputar.globalsub;.section/pmprogram;sub: se=-1; sr=ashiftax0(hi); ax0=sr1;ax1=sr0; sr=ashiftay0(hi); ay0=sr1;ay1=sr0; ar=ax1-ay1; ar=ax0-ay0+c-1; rts;//filename:multi.asm//purpose:multiplyopcode.globalmulti;.section/pmprogram;//(a+bi)*(c+di)=ac-bd+(ad+bc)imulti: mr=mx0*my0(ss); mr=mr-mx1*my1(ss); sr=mx0*my1(ss); sr=sr+mx1*my0(ss); rts;//filename:halfbutterfly.asm//purpose:halfbutterflyopcode.globalhalfbutterfly;.section/pmprogram;halfbutterfly: cntr=4;//doitfor4timesdocircleuntilce; ax0=dm(i0+=m0);//readIdata ax1=dm(i1+=m0);//readQdata my0=dm(i2+=m3);//readcosdata my1=dm(i3+=m3);//readsindata mx0=ax0;//formultiply mx1=ax1;//formultiply callmulti; mx0=mr1;//Idata mx1=sr1;//Qdata dm(i4+m5)=mx0;//saveIdata dm(i5+m5)=mx1;//saveQdata se=-1; sr=ashiftax0(hi); ax0=sr1; dm(i4+=m4)=ax0;//saveIdata sr=ashiftax1(hi); ax1=sr1;circle: dm(i5+=m4)=ax1;//saveQdata rts;//filename:totalbutterfly.asm//purpose:totalbutterflyopcode.globaltotalbutterfly;.section/pmprogram;totalbutterfly:cntr=2;//do2timesdocircle2untilce; m2=0; my0=dm(i2+=m3);//readcosdata my1=dm(i3+=m3);//readsindata //ar=ax0-ay0;//realI-realQ ay0=dm(i0+m1);//readIdata ax0=dm(i0+m2);//readIdata callsub; mx0=ar; //ar=ax1-ay1;//imageI-imageQ ay0=dm(i1+m1);//readQdata ax0=dm(i1+m2);//readQdata callsub; mx1=ar; callmulti; dm(i4+m5)=mr1;//saveIdata dm(i5+m5)=sr1;//saveQdata //calladd; //ar=ax0+ay0;//realI+realQ ay0=dm(i0+m1);//readIdata ax0=dm(i0+=m0);//readIdata calladd; se=-1; s