程序员专用DSP算法：数学原理与源码实现（(Dsp)AlgorithmsForProgrammers）摘要：《(Dsp)AlgorithmsForProgrammers》是一本面向程序员的DSP算法与数学原理实战手册，聚焦数字信号处理核心算法的数学原理与代码实现。全书系统讲解离散傅里叶变换、基2/基4/分裂基FFT、卷积、Z变换、哈特利变换、数论变换、沃尔什变换、哈尔变换等核心算法，给出递归与非递归实现、三角函数计算优化、内存访问改进及代码自动生成方法。同时涵盖位操作、排列、排序、组合算法、高精度运算等底层技术，提供可直接复用的伪代码与C++源码，兼顾理论推导与工程效率，适合DSP开发、算法实现与高性能计算参考。AlgorithmsforprogrammersideasandsourcecodeThisdocumentisworkinprogress:readthe”importantremarks”nearthebeginningJ¨orgArndtarndt@jjj.deThisdocument1wasLATEX’datSeptember26,2002Someimportantremarksaboutthisdocument.Thisdraftisintendedtoturnintoabookaboutselectedalgorithms.Theaudienceinmindarepro-grammerswhoareinterestedinthetreatedalgorithmsandactuallywanttohave/createworkingandreasonablyoptimizedcode.Theprintablefullversionwillalwaysstayonlineforfreedownload.ItisplannedtoalsomakepartsoftheTEXsources(plusthescriptsusedforautomation)available.RightnowafewfilesoftheTEXsourcesandallextractedpseudo-codesnippets1areonline.TheC++-sourcesareonlineaspartofFXTorhfloat(arithmeticalalgorithms).ThequalityandspeedofdevelopmentdoesdependonthefeedbackthatIreceivefromyou.Yourcriticismconcerninglanguage,style,correctness,omissions,technicalitiesandeventhegoalssethereisverywelcome.Thankstothose2whohelpedtoimprovethisdocumentsofar!Thanksalsotothepeoplewhosharetheirideas(orsourcecode)onthenet.Itrytogiveduereferencestooriginalsources/authorswhereverIcan.However,IaminnowayanexpertforhistoryofalgorithmsandIprettysurewillneverbeone.Soifyoufeelthatareferenceismissingsomewhere,letmeknow.Newchapters/sectionsappearassoonastheycontainanythinguseful,sometimesjustlistingsorremarksoutliningwhatistoappearthere.A”TBD:somethingtobedone”isaremindertomyselftofillinsomethingthatismissingorwouldbenicetohave.ThestylevariesfromchaptertochapterwhichIdonotconsiderbadperse:whilesometopics(e.g.FFTs)needaclearandexplicitintroductionothers(e.g.thebitwizardrychapter)seemtobebestpresentedbybasicallyshowingthecodewithjustafewcomments.Stillotherparts(e.g.sorting)arepresentedelsewhereextremelywellsoIwillintroducethebasicideasonlyveryshortlyandsupplysome(hopefully)usefulcode.Sprachewillpartlygoaway:using/includingtheactualcodefromFXTwillbebeneficialtoboththisdocumentandFXTitself.Thegoalistoautomaticallyincludethefunctionsreferenced.Clearly,thiswilldrasticallyreducethechanceoferrorsintheshowncode(andatthesametimedrasticallyreducetheworkloadforme).InitiallyIplannedtowriteaninterpreterforSprache,itjustneverhappened.AtthesametimeFXTwillbebetterdocumentedwhichitreallyneeds.AsaconsequenceSprachewillonlybeusedwhenthereisaclearadvantagetodoso,mainlywhenthecorrespondingC++doesnotappeartobeselfexplanatory.LargerpiecesofcodewillbepresentedinC++.AtinystarteraboutC++(somegoodreasonsinfavorofC++andsomeoftheverybasicsofclasses/overloading/templates)willbeincluded.Cprogrammersdonotneedtobeshockedbytheζ++’:onlyanratherminimalsetoftheC++featuresisused.Thetheorem-likeenvironmentforthecodesshallcompletelygoaway.Itleadstoduplicationofstate-ments,especiallywithnon-pseudocode(runningtext,descriptionintheenvironmentandcommentsatthebeginoftheactualcode).Enjoyreading!ListofimportantSymbols沦xrealpartofx$ximaginarypartofxx∗complexconjugateofxaasequence,e.g.{a0,a1,...,an−1},theindexalwaysstartswithzero.atransformed(e.g.Fouriertransformed)sequencememphasizethatthesequencestotheleftandrightarealloflengthm=F[a](=c)(discrete)Fouriertransform(FT)ofa,ck=Σx(n)0(1)axzxkwherez=e±2πi/nF−1[a]inverse(discrete)Fouriertransform(IFT)ofa,F−1[a]k=Σx(n)0(1)axz−xkSkaasequencecwithelementscx:=axe±k2πix/nH[a]discreteHartleytransform(HT)ofaasequencereversedaroundelementwithindexn/2aSthesymmetricpartofasequence:aS:=a+aaAtheantisymmetricpartofasequence:aA:=a−aZ[a]discretez-transform(ZT)ofaWv[a]discreteweightedtransformofa,weight(sequence)vW1[a]inversediscreteweightedtransformofa,weightva④bcyclic(orcircular)convolutionofsequenceawithsequenceba④acbacyclic(orlinear)convolutionofsequenceawithsequenceba④−bnegacyclic(orskewcircular)convolutionofsequenceawithsequenceba④{v}bweightedconvolutionofsequenceawithsequenceb,weightva④⊕bdyadicconvolutionofsequenceawithsequencebn\NndividesNn⊥mgcd(n,m)=1a(j%m)sequenceconsistingoftheelementsofawithindicesk:k≡jmodme.g.a(even),a(odd)a(0%2),a(1%2)a(j/m)sequenceconsistingoftheelementsofawithindicesk:j·n/m≤k<(j+1)·n/me.g.a(left),a(right)a(0/2),a(1/2)Chapter1TheFouriertransform1.1ThediscreteFouriertransformThediscreteFouriertransform(DFTorsimplyFT)ofacomplexsequenceaoflengthnisdefinedaszisann-throotofunity:zn=1.Backtransform(orinversediscreteFouriertransformIDFTorsimplyIFT)isthenToseethis,considerelementyoftheIFToftheFTofa:AsΣk(z①__y)k=nforx=yandzeroelse(becausezisann-throotofunity).ThereforethewholeexpressionisequaltowhereHerewewillcalltheFTwiththeplusintheexponenttheforwardtransform.Thechoiceisactuallyarbitrary1.CHAPTER1.THEFOURIERTRANSFORMThenormalizationfactorinfrontoftheFTsumsissometimesreplacedbyasingleinfrontoftheinverseFTsumwhichisoftenconvenientincomputation.Then,ofcourse,Parseval’sequationhastobemodifiedaccordingly.AstraightforwardimplementationofthediscreteFouriertransform,i.e.thecomputationofnsumseachoflengthnrequires~n2operations:[FXT:slowftinslow/slowft.cc]ismustbe+1(forwardtransform)or__1(backwardtransform),SinCos(x)returnsaComplex(cos(x),sin(x)).AfastFouriertransform(FFT)algorithmisanalgorithmthatimprovestheoperationcounttopropor-tionalnΣk__1),wheren=p1p2···pmisafactorizationofn.Incaseofapowern=pmthevaluecomputeston(p__1)logp(n).Inthespecialcasep=2evenn/2log2(n)(complex)multiplicationssuffice.ThereareseveraldifferentFFTalgorithmswithmanyvariants.1.2SymmetriesoftheFouriertransformAbitofnotationturnsouttobeuseful:Letabethesequencea(lengthn)reversedaroundelementwithindexn/2:LetaS,aAbethesymmetric,antisymmetricpartofthesequencea,respectively(Theelementswithindices0andn/2ofaAarezero).Nowleta∈R(meaningthateachelementofais∈R),thenTheFTofarealsymmetricsequenceisrealandsymmetricandtheFTofarealantisymmetricsequenceispurelyimaginaryandantisymmetric.TherebytheFTofageneralrealsequenceisthecomplexconjugateofitsreversed:Similarly,forapurelyimaginarysequenceb∈iR:TheFTofacomplexsymmetric/antisymmetricsequenceissymmetric/antisymmetric,respectively.1.3Radix2FFTalgorithms1.3.1AlittlebitofnotationAlwaysassumeaisalength-nsequence(napoweroftwo)inwhatfollows:Leta(even),a(odd)denotethe(length-n/2)subsequencesofthoseelementsofathathaveevenoroddindices,respectively.Leta(left)denotethesubsequenceofthoseelementsofathathaveindices0...n/2__1.Similarly,a(right)forindicesn/2n__1.LetSkadenotethesequencewithelementsaxe§k2πix/nwherenisthelengthofthesequenceaandthesignisthatofthetransform.ThesymbolSshallsuggestashiftoperator.InthenexttwosectionsonlyS1/2willappear.S0istheidentityoperator.1.3.2Decimationintime(DIT)FFTThefollowingobservationisthekeytothedecimationintime(DIT)FFT2algorithm:Forneventhek-thelementoftheFouriertransformiswherez=e§i2π/nandk∈{0,1,,n__1}.Thelastidentitytellsushowtocomputethek-thelementofthelength-nFouriertransformfromthelength-n/2Fouriertransformsoftheevenandoddindexedsubsequences.Toactuallyrewritethelength-nFTintermsoflength-n/2FTsonehastodistinguishthecases0≤k<n/2andn/2≤k<n,thereforewerewritek∈{0,1,2,...,n__1}ask=j+δwherej∈{0,1,...,n/2__1},δ∈{0,1}.Notingthatz2isjusttherootofunitythatappearsinalength-n/2FTonecanrewritethelasttwoequationsastheIdea1.1(FFTradix2DITstep)Radix2decimationintimestepfortheFFT:(Hereitissilentlyassumedthat’+’or’__’betweentwosequencesdenoteselementwiseadditionorsubtraction.)Thelength-ntransformhasbeenreplacedbytwotransformsoflengthn/2.Ifnisapowerof2thisschemecanbeappliedrecursivelyuntillength-onetransforms(identityoperation)arereached.Therebytheoperationcountisimprovedtoproportionaln·log2(n):Therearelog2(n)splittingsteps,theworkineachstepisproportionalton.Code1.1(recursiveradix2DITFFT)Pseudocodeforarecursiveprocedureofthe(radix2)DITFFTalgorithm,ismustbe+1(forwardtransform)or-1(backwardtransform):Thedatalengthnmustbeapowerof2.Theresultisinx[].Notethatnormalization(i.e.multiplicationofeachelementofx[]by1/√n)isnotincludedhere.[FXT:recursivedit2fftinslow/recfft2.cc]TheprocedureusesthesubroutineCode1.2(Fouriershift)Foreachelementinc[0..n-1]replacec[k]byc[k]timesev2πik/n.UsedwithU=§1/2fortheFouriertransform.TherecursiveFFT-procedureinvolvesnlog2(n)functioncalls,whichcanbeavoidedbyrewritingitinanon-recursiveway.Onecanevendoalloperationsinplace,notemporaryworkspaceisneededatall.Thepriceisthenecessityofanadditionaldatareordering:Theprocedurerevbin_permute(a[],n)rearrangesthearraya[]inawaythateachelementaxisswappedwitha˜x,where˜xisobtainedfromxbyreversingitsbinarydigits.Thisisdiscussedinsection8.1.Code1.3(radix2DITFFT,localized)Pseudocodeforanon-recursiveprocedureofthe(radix2)DITalgorithm,ismustbe-1or+1:[FXT:dit2-fft-localizedinfft/fftdit2.cc]Thisversionofanon-recursiveFFTprocedurealreadyavoidsthecallingoverheadanditworksinplace.Itworksasgiven,butisabitwasteful.The(expensive!)computatione:=exp(is*2*PI*I*j/m)isdonen/2·log2(n)times.Toreducethenumberoftrigonometriccomputations,onecansimplyswapthetwoinnerloops,leadingtothefirstζrealworld’FFTprocedurepresentedhere:Code1.4(radix2DITFFT)Pseudocodeforanon-recursiveprocedureofthe(radix2)DITalgo-rithm,ismustbe-1or+1:procedurefft_dit2(a[],ldn,is)//complexa[0..2**ldn-1]input,resultSwappingthetwoinnerloopsreducesthenumberoftrigonometric(exp())computationstonbutleadstoafeaturethatmanyFFTimplementationsshare:Memoryaccessishighlynonlocal.Foreachrecursionstage(valueofldm)thearrayistraversedmhtimeswithn/maccessesinstridesofmh.Asmhisapowerof2thiscan(oncomputersthatusememorycache)haveaverynegativeperformanceimpactforlargevaluesofn.OnacomputerwheretheCPUclock(366MHz,AMDK6/2)is5.5timesfasterthanthememoryclock(66MHz,EDO-RAM)IfoundthatindeedforsmallnthelocalizedFFTisslowerbyafactorofabout0.66,butforlargenthesameratioisinfavouroftheζnaive’procedure!Itisagoodideatoextracttheldm==1stageoftheoutermostloop,thisavoidscomplexmultiplicationswiththetrivialfactors1+0i:Replace更多内容见附件免责声明：本页面/内容部分素材来源于互联网公开信息，旨在传递更多信息，不代表本平台立场。版权归原作者或机构所有，如涉及侵权，请通过平台联系我们，我们将在核实后第一时间处理。本平台对转载内容的真实性、准确性不作任何保证，用户需自行判断并承担使用风险
