Goertzel algorithm好文收集
2008-08-02 11:45
246 查看
Goertzelalgorithm
FromWikipedia,thefreeencyclopedia
Jumpto:navigation,search
TheGoertzelalgorithmisadigitalsignalprocessing(DSP)techniqueforidentifyingfrequencycomponentsofasignal,publishedbyDr.GeraldGoertzelin1958.WhilethegeneralFastFouriertransform(FFT)algorithmcomputesevenlyacrossthebandwidthoftheincomingsignal,theGoertzelalgorithmlooksatspecific,predeterminedfrequency.
Apracticalapplicationofthisalgorithmisthatofrecognizingthetonesproducedbythebuttonspushedonatelephonekeypad.ThisapplicationisillustratedbyanimplementationofthealgorithminCasshownbelowtoproduceaDTMFtonedetector.
//if(window.showTocToggle){vartocShowText="show";vartocHideText="hide";showTocToggle();}
//]]>
[edit]Explanationofalgorithm
TheGoertzelalgorithmcomputesasequence,s(n),givenaninputsequence,x(n),as
s(n)=x(n)+2cos(2πω)s(n−1)−s(n−2)
wheres(−2)=s(−1)=0andωissomefrequencyofinterest,incyclespersample,whichshouldbelessthan1/2.Thiseffectivelyimplementsasecond-orderIIRfilterwithpolesate+2πiωande−2πiω,andrequiresonlyonemultiplication(assuming2cos(2πω)ispre-computed),oneadditionandonesubtractionperinputsample.Forrealinputs,theseoperationsarereal.
TheZtransformofthisprocessis
Applyinganadditional,FIR,transformoftheform
willgiveanoveralltransformof
Thetime-domainequivalentofthisoveralltransformis
which,whenx(n)=0foralln<0,becomes
or,theequationforthe(n+1)-sampleDFTofx,evaluatedforωandmultipliedbythescalefactore+2πiωn.
NoticethatapplyingtheadditionaltransformY(z)/S(z)onlyrequiresthelasttwosamplesofthessequence.Consequently,uponprocessingNsamplesx(0)...x(N−1),thelasttwosamplesfromthessequencecanbeusedtocomputethevalueofaDFTbinwhichcorrespondstothechosenfrequencyωas
X(ω)=y(N−1)e−2πiω(N−1)=(s(N−1)−e−2πiωs(N−2))e−2πiω(N−1)
ForthespecialcaseoftenfoundwhencomputingDFTbins,whereωN=kforsomeinteger,k,thissimplifiesto
X(ω)=(s(N−1)−e−2πiωs(N−2))e+2πiω=e+2πiωs(N−1)−s(N−2)
Ineithercase,thecorrespondingpowercanbecomputedusingthesamecosinetermrequiredtocomputesas
X(ω)X'(ω)=s(N−2)2+s(N−1)2−2cos(2πω)s(N−2)s(N−1)
Whenimplementedinageneral-purposeprocessor,valuesfors(n−1)ands(n−2)canberetainedinvariablesandnewvaluesofscanbeshiftedthroughastheyarecomputed,assumingthatonlythefinaltwovaluesofthessequencearerequired.Thecodemaythenbeasfollows:
s_prev=0
s_prev2=0
coeff=2*cos(2*PI*normalized_frequency);
foreachsample,x
,
s=x
+coeff*s_prev-s_prev2;
s_prev2=s_prev;
s_prev=s;
end
power=s_prev2*s_prev2+s_prev*s_prev-coeff*s_prev2*s_prev;
[edit]Computationalcomplexity
InordertocomputeasingleDFTbinforacomplexsequenceoflengthN,thisalgorithmrequires2Nmultipliesand4Nadd/subtractoperationswithintheloop,aswellas4multipliesand4add/subtractoperationstocomputeX(ω),foratotalof2N+4multipliesand4N+4add/subtractoperations(forrealsequences,therequiredoperationsarehalfthatamount).Incontrast,theFastFouriertransform(FFT)requires2log2Nmultipliesand3log2Nadd/subtractoperationsperDFTbin,butmustcomputeallNbinssimultaneously(similaroptimizationsareavailabletohalvethenumberofoperationsinanFFTwhentheinputsequenceisreal).
WhenthenumberofdesiredDFTbins,M,issmall(e.g.,whendetectingDTMFtones),itiscomputationallyadvantageoustoimplementtheGoertzelalgorithm,ratherthantheFFT.Approximately,thisoccurswhen
orif,forsomereason,Nisnotanintegralpowerof2whileyousticktoasimpleFFTalgorithmlikethe2-radixCooley-TukeyFFTalgorithm,andzero-paddingthesamplesouttoanintegralpowerof2wouldviolate
Moreover,theGoertzelalgorithmcanbecomputedassamplescomein,andtheFFTalgorithmmayrequirealargetableofNpre-computedsinesandcosinesinordertobeefficient.
Ifmultiplicationsarenotconsideredasdifficultasadditions,orviceversa,the5/6ratiocanshiftbetweenanythingfrom3/4(additionsdominate)to1/1(multiplicationsdominate).
[edit]Practicalconsiderations
ThetermDTMForDual-ToneMultiFrequencyistheofficialnameofthetonesgeneratedfromatelephonekeypad.(AT&Tusedthetrademark"Touch-Tone"foritsDTMFdialingservice.[1])TheoriginalkeypadsweremechanicalswitchestriggeringRCcontrolledoscillators.[citationneeded ]Thedigitdetectorswerealsotunedcircuits.TheinterestindecodingDTMFishighbecauseofthelargenumbersofphonesgeneratingthesetypesoftones.
Atpresent,DTMFdetectorsaremostoftenimplementedasnumericalalgorithmsoneithergeneralpurposecomputersoronfastdigitalsignalprocessors.Thealgorithmshownbelowisanexampleofsuchadetector.
However,thisalgorithmneedsanadditionalpost-processingsteptocompletelyimplementafunctionalDTMFtonedetector.DTMFtoneburstscanbeasshortas50milli-secondsoraslongasseveralseconds.Thetoneburstcanhavenoiseordropoutswithinitwhichmustbeignored.TheGoertzelalgorithmproducesmultipleoutputs;apost-processingstepneedstosmooththeseoutputsintooneoutputpertoneburst.
Oneadditionalproblemisthatthealgorithmwillsometimesproducespuriousoutputsbecauseofawindowperiodthatisnotcompletelyfilledwithsamples.ImagineaDTMFtoneburstandthenimaginethewindowsuperimposedoverthistoneburst.Obviously,thedetectorisrunningatafixedrateandthetoneburstisnotguaranteedtoarrivealignedwiththetimingofthedetector.Sosomewindowintervalsontheleadingandtrailingedgesofthetoneburstwillnotbeentirelyfilledwithvalidtonesamples.Worse,RC-basedtonegeneratorswilloftenhavevoltagesag/surgerelatedanomaliesattheleadingandtrailingedgesofthetoneburst.Thesealsocancontributetospuriousoutputs.
Itishighlylikelythatthisdetectorwillreportfalseorincorrectresultsattheleadingandtrailingedgesofthetoneburstduetoalackofsufficientvalidsampleswithinthewindow.Inaddition,thetonedetectormustbeabletotoleratetonedropoutswithinthetoneburstandthesecanproduceadditionalfalsereportsduetothesamewindowingeffects.
Thepost-processingsystemcanbeimplementedasastatisticalaggregatorwhichwillexamineaseriesofoutputsofthealgorithmbelow.Thereshouldbeacounterforeachpossibleoutput.Theseallstartoutatzero.Thedetectorstartsproducingoutputsanddependingontheoutput,theappropriatecounterisincremented.Finally,thedetectorstopsgeneratingoutputsforlongenoughthatthetoneburstcanbeconsideredtobeover.ThecounterwiththehighestvaluewinsandshouldbeconsideredtobetheDTMFdigitsignaledbythetoneburst.
WhileitistruethatthereareeightpossiblefrequenciesinaDTMFtone,thealgorithmasoriginallyenteredonthispagewascomputingafewmorefrequenciessoastohelprejectfalsetones(talkoff).Noticethepeaktonecounterloop.Thischeckstoseethatonlytwotonesareactive.Ifmorethanthisarefoundthenthetoneisrejected.
[edit]SamplecodeforaDTMFdetector
#defineSAMPLING_RATE8000
#defineMAX_BINS8
#defineGOERTZEL_N92
intsample_count;
doubleq1[MAX_BINS];
doubleq2[MAX_BINS];
doubler[MAX_BINS];
/*
*coef=2.0*cos((2.0*PI*k)/(float)GOERTZEL_N));
*Wherek=(int)(0.5+((float)GOERTZEL_N*target_freq)/SAMPLING_RATE));
*
*Moresimply:coef=2.0*cos((2.0*PI*target_freq)/SAMPLING_RATE);
*/
doublefreqs[MAX_BINS]=
{
697,
770,
852,
941,
1209,
1336,
1477,
1633
};
doublecoefs[MAX_BINS];
/*----------------------------------------------------------------------------
*calc_coeffs
*----------------------------------------------------------------------------
*Thisiswherewecalculatethecorrectco-efficients.
*/
voidcalc_coeffs()
{
intn;
for(n=0;n<MAX_BINS;n++)
{
coefs[n]=2.0*cos(2.0*3.141592654*freqs[n]/SAMPLING_RATE);
}
}
/*----------------------------------------------------------------------------
*post_testing
*----------------------------------------------------------------------------
*Thisiswherewelookatthebinsanddecideifwehaveavalidsignal.
*/
voidpost_testing()
{
introw,col,see_digit;
intpeak_count,max_index;
doublemaxval,t;
inti;
char*row_col_ascii_codes[4][4]={
{"1","2","3","A"},
{"4","5","6","B"},
{"7","8","9","C"},
{"*","0","#","D"}};
/*Findthelargestintherowgroup.*/
row=0;
maxval=0.0;
for(i=0;i<4;i++)
{
if(r>maxval)
{
maxval=r[i];
row=i;
}
}
[i]/*Findthelargestinthecolumngroup.*/
col=4;
maxval=0.0;
for(i=4;i<8;i++)
{
if(r>maxval)
{
maxval=r[i];
col=i;
}
}
[i]/*Checkforminimumenergy*/
if(r[row]<4.0e5)/*2.0e5...1.0e8nochange*/
{
/*energynothighenough*/
}
elseif(r[col]<4.0e5)
{
/*energynothighenough*/
}
else
{
see_digit=TRUE;
/*Twistcheck
*CEPT=>twist<6dB
*AT&T=>forwardtwist<4dBandreversetwist<8dB
*-ndB<10log10(v1/v2),wherev1<v2
*-4dB<10log10(v1/v2)
*-0.4<log10(v1/v2)
*0.398<v1/v2
*0.398*v2<v1
*/
if(r[col]>r[row])
{
/*Normaltwist*/
max_index=col;
if(r[row]<(r[col]*0.398))/*twist>4dB,error*/
see_digit=FALSE;
}
else/*if(r[row]>r[col])*/
{
/*Reversetwist*/
max_index=row;
if(r[col]<(r[row]*0.158))/*twist>8db,error*/
see_digit=FALSE;
}
/*Signaltonoisetest
*AT&Tstatesthatthenoisemustbe16dBdownfromthesignal.
*Herewecountthenumberofsignalsabovethethresholdand
*thereoughttobeonlytwo.
*/
if(r[max_index]>1.0e9)
t=r[max_index]*0.158;
else
t=r[max_index]*0.010;
peak_count=0;
for(i=0;i<8;i++)
{
if(r>t)
peak_count++;
}
if(peak_count>2)
see_digit=FALSE;
if(see_digit)
{
printf("%s",row_col_ascii_codes[row][col-4]);
fflush(stdout);
}
}
}
/*----------------------------------------------------------------------------[i]
*goertzel
*----------------------------------------------------------------------------
*/
voidgoertzel(intsample)
{
doubleq0;
ui32i;
sample_count++;
/*q1[0]=q2[0]=0;*/
for(i=0;i<MAX_BINS;i++)
{
q0=coefs*q1[i]-q2[i]+sample;
q2[i]=q1[i];
q1[i]=q0;
}
if(sample_count==GOERTZEL_N)
{
for(i=0;i<MAX_BINS;i++)
{
r[i]=(q1[i]*q1[i])+(q2[i]*q2[i])-(coefs[i]*q1[i]*q2[i]);
q1[i]=0.0;
q2[i]=0.0;
}
post_testing();
sample_count=0;
}
}
[edit]References
^ USPTOtrademarkentryforTouch-Tone,retrievedonMarch29,2007.
[edit]Externallinks
http://ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html
http://www.embedded.com/story/OEG20020819S0057
http://www.embedded.com/showArticle.jhtml?articleID=17301593
http://www.numerix-dsp.com/goertzel.html
http://www.mathworks.com/access/helpdesk/help/toolbox/signal/goertzel.html
DTMFdecodingwitha1-bitA/Dconverter
ModifiedGoertzelAlgorithminDTMFDetectionUsingtheTMS320C80DSP
TheGoertzelAlgorithm
TheGoertzelalgorithmcanperformtonedetectionusingmuchlessCPUhorsepowerthantheFastFourierTransform,butmanyengineershaveneverheardofit.Thisarticleattemptstochangethat.
MostengineersarefamiliarwiththeFastFourierTransform(FFT)andwouldhavelittletroubleusinga"canned"FFTroutinetodetectoneormoretonesinanaudiosignal.Whatmanydon'tknow,however,isthatifyouonlyneedtodetectafewfrequencies,amuchfastermethodisavailable.It'scalledtheGoertzelalgorithm.
Tonedetection
Manyapplicationsrequiretonedetection,suchas:
DTMF(touchtone)decoding
Callprogress(dialtone,busy,andsoon)decoding
Frequencyresponsemeasurements(sendingatonewhilesimultaneouslyreadingbacktheresult)-ifyoudothisforarangeoffrequencies,theresultingfrequencyresponsecurvecanbeinformative.Forexample,thefrequencyresponsecurveofatelephonelinetellsyouifanyloadcoils(inductors)arepresentonthatline.
AlthoughdedicatedICsexistfortheapplicationsabove,implementingthesefunctionsinsoftwarecostsless.Unfortunately,manyembeddedsystemsdon'thavethehorsepowertoperformcontinuousreal-timeFFTs.That'swheretheGoertzelalgorithmcomesin.
Inthisarticle,IdescribewhatIcallabasicGoertzelandanoptimizedGoertzel.
ThebasicGoertzelgivesyourealandimaginaryfrequencycomponentsasaregularDiscreteFourierTransform(DFT)orFFTwould.Ifyouneedthem,magnitudeandphasecanthenbecomputedfromthereal/imaginarypair.
TheoptimizedGoertzelisevenfaster(andsimpler)thanthebasicGoertzel,butdoesn'tgiveyoutherealandimaginaryfrequencycomponents.Instead,itgivesyoutherelativemagnitudesquared.Youcantakethesquarerootofthisresulttogettherelativemagnitude(ifneeded),butthere'snowaytoobtainthephase.
Inthisshortarticle,Iwon'ttrytoexplainthetheoreticalbackgroundofthealgorithm.Idogivesomelinksattheendwhereyoucanfindmoredetailedexplanations.Icantellyouthatthealgorithmworkswell,havinguseditinallofthetonedetectionapplicationspreviouslylisted(andothers).
AbasicGoertzel
Firstaquickoverviewofthealgorithm:someintermediateprocessingisdoneineverysample.TheactualtonedetectionoccurseveryNthsample.(I'lltalkmoreaboutNinaminute.)
AswiththeFFT,youworkwithblocksofsamples.However,thatdoesn'tmeanyouhavetoprocessthedatainblocks.Thenumericalprocessingisshortenoughtobedoneintheveryinterruptserviceroutine(ISR)thatisgatheringthesamples(ifyou'regettinganinterruptpersample).Or,ifyou'regettingbuffersofsamples,youcangoaheadandprocessthemabatchatatime.
BeforeyoucandotheactualGoertzel,youmustdosomepreliminarycalculations:
Decideonthesamplingrate.
Choosetheblocksize,N.
Precomputeonecosineandonesineterm.
Precomputeonecoefficient.
Thesecanallbeprecomputedonceandthenhardcodedinyourprogram,savingRAMandROMspace;oryoucancomputethemon-the-fly.
Samplingrate
Yoursamplingratemayalreadybedeterminedbytheapplication.Forexample,intelecomapplications,it'scommontouseasamplingrateof8kHz(8,000samplespersecond).Alternatively,youranalog-to-digitalconverter(orCODEC)mayberunningfromanexternalclockorcrystaloverwhichyouhavenocontrol.
Ifyoucanchoosethesamplingrate,theusualNyquistrulesapply:thesamplingratewillhavetobeatleasttwiceyourhighestfrequencyofinterest.Isay"atleast"becauseifyouaredetectingmultiplefrequencies,it'spossiblethatanevenhighersamplingfrequencywillgivebetterresults.Whatyoureallywantisforeveryfrequencyofinteresttobeanintegerfactorofthesamplingrate.
Blocksize
GoertzelblocksizeNislikethenumberofpointsinanequivalentFFT.Itcontrolsthefrequencyresolution(alsocalledbinwidth).Forexample,ifyoursamplingrateis8kHzandNis100samples,thenyourbinwidthis80Hz.
ThiswouldsteeryoutowardsmakingNashighaspossible,togetthehighestfrequencyresolution.ThecatchisthatthehigherNgets,thelongerittakestodetecteachtone,simplybecauseyouhavetowaitlongerforallthesamplestocomein.Forexample,at8kHzsampling,itwilltake100msfor800samplestobeaccumulated.Ifyou'retryingtodetecttonesofshortduration,youwillhavetousecompatiblevaluesofN.
ThethirdfactorinfluencingyourchoiceofNistherelationshipbetweenthesamplingrateandthetargetfrequencies.Ideallyyouwantthefrequenciestobecenteredintheirrespectivebins.Inotherwords,youwantthetargetfrequenciestobeintegermultiplesofsample_rate/N.
Thegoodnewsisthat,unliketheFFT,Ndoesn'thavetobeapoweroftwo.
Precomputedconstants
Onceyou'veselectedyoursamplingrateandblocksize,it'sasimplefive-stepprocesstocomputetheconstantsyou'llneedduringprocessing:
[i]w =(2*π/N)*k
cosine=cosw
sine=sinw
coeff=2*cosine
Fortheper-sampleprocessingyou'regoingtoneedthreevariables.Let'scallthemQ0,Q1,andQ2.
Q1isjustthevalueofQ0lasttime.Q2isjustthevalueofQ0twotimesago(orQ1lasttime).
Q1andQ2mustbeinitializedtozeroatthebeginningofeachblockofsamples.Foreverysample,youneedtorunthefollowingthreeequations:
Q0=coeff*Q1-Q2+sample
Q2=Q1
Q1=Q0
Afterrunningtheper-sampleequationsNtimes,it'stimetoseeifthetoneispresentornot.
real=(Q1-Q2*cosine)
imag=(Q2*sine)
magnitude2=real2+imag2
Asimplethresholdtestofthemagnitudewilltellyouifthetonewaspresentornot.ResetQ2andQ1tozeroandstartthenextblock.
AnoptimizedGoertzel
TheoptimizedGoertzelrequireslesscomputationthanthebasicone,attheexpenseofphaseinformation.
Theper-sampleprocessingisthesame,buttheendofblockprocessingisdifferent.Insteadofcomputingrealandimaginarycomponents,andthenconvertingthoseintorelativemagnitudesquared,youdirectlycomputethefollowing:
magnitude2=Q12+Q22-Q1*Q2*coeff
ThisistheformofGoertzelI'veusedmostoften,anditwasthefirstoneIlearnedabout.
Pullingitalltogether
Listing1showsashortdemoprogramIwrotetoenableyoutotest-drivethealgorithm.ThecodewaswrittenandtestedusingthefreelyavailableDJGPPC/C++compiler.Youcanmodifythe#definesnearthetopofthefiletotryoutdifferentvaluesofN,sampling_rate,andtarget_frequency.
Theprogramdoestwodemonstrations.Inthefirstone,bothformsoftheGoertzelalgorithmareusedtocomputerelativemagnitudesquaredandrelativemagnitudeforthreedifferentsynthesizedsignals:onebelowthetarget_frequency,oneequaltothetarget_frequency,andoneabovethetarget_frequency.
You'llnoticethattheresultsarenearlyidentical,regardlessofwhichformoftheGoertzelalgorithmisused.You'llalsonoticethatthedetectorvaluespeaknearthetargetfrequency.
Intheseconddemonstration,asimulatedfrequencysweepisrun,andtheresultsofjustthebasicGoertzelareshown.Again,you'llnoticeaclearpeakinthedetectoroutputnearthetargetfrequency.Figure1showstheoutputofthecodeinListing1.
Figure1:Demooutput
ForSAMPLING_RATE=8000.000000N=205andFREQUENCY=941.000000,k=24andcoeff=1.482867
Toavoidfalsedetectionsinproductioncode,youwillprobablywanttoqualifytherawdetectorreadingsbyusingadebouncingtechnique,suchasrequiringmultipledetectionsinarowbeforereportingatone'spresencetotheuser.
Asyoucansee,theGoertzelalgorithmdeservestobeaddedtoyoursignalprocessingtoolbox.
KevinBankshasbeendevelopingembeddedsystemsfor19years,asaconsultantandasanemployeeatcompaniesincludingSCI,TxPORT,DiscoveryCom,andNokia.Currentlyheisbackinconsultingmode.Hecanbereachedatkbanks@hiwaay.net.
References
Herearesomelinkstootherresourcesyoumayfinduseful:
www.numerix-dsp.com/goertzel.html
ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html
www.analogdevices.com/library/dspManuals/Using_ADSP-2100_Vol1_books.html
www.analogdevices.com/library/dspManuals/pdf/Volume1/Chapter_14.pdf
www.analogdevices.com/library/dspManuals/Using_ADSP-2100_Vol2_books.html
www.analogdevices.com/library/dspManuals/pdf/2100Chapter_8.pdf
www.delorie.com/djgpp/
FromWikipedia,thefreeencyclopedia
Jumpto:
TheGoertzelalgorithmisa
Apracticalapplicationofthisalgorithmisthatofrecognizingthetonesproducedbythebuttonspushedonatelephonekeypad.ThisapplicationisillustratedbyanimplementationofthealgorithminCasshownbelowtoproducea
| Contents [ |
//if(window.showTocToggle){vartocShowText="show";vartocHideText="hide";showTocToggle();}
//]]>
[
TheGoertzelalgorithmcomputesasequence,s(n),givenaninputsequence,x(n),as
s(n)=x(n)+2cos(2πω)s(n−1)−s(n−2)
wheres(−2)=s(−1)=0andωissomefrequencyofinterest,incyclespersample,whichshouldbelessthan1/2.Thiseffectivelyimplementsasecond-order
The
Applyinganadditional,
willgiveanoveralltransformof
Thetime-domainequivalentofthisoveralltransformis
which,whenx(n)=0foralln<0,becomes
or,theequationforthe(n+1)-sample
NoticethatapplyingtheadditionaltransformY(z)/S(z)onlyrequiresthelasttwosamplesofthessequence.Consequently,uponprocessingNsamplesx(0)...x(N−1),thelasttwosamplesfromthessequencecanbeusedtocomputethevalueofa
X(ω)=y(N−1)e−2πiω(N−1)=(s(N−1)−e−2πiωs(N−2))e−2πiω(N−1)
Forthespecialcaseoftenfoundwhencomputing
X(ω)=(s(N−1)−e−2πiωs(N−2))e+2πiω=e+2πiωs(N−1)−s(N−2)
Ineithercase,thecorrespondingpowercanbecomputedusingthesamecosinetermrequiredtocomputesas
X(ω)X'(ω)=s(N−2)2+s(N−1)2−2cos(2πω)s(N−2)s(N−1)
Whenimplementedinageneral-purposeprocessor,valuesfors(n−1)ands(n−2)canberetainedinvariablesandnewvaluesofscanbeshiftedthroughastheyarecomputed,assumingthatonlythefinaltwovaluesofthessequencearerequired.Thecodemaythenbeasfollows:
s_prev=0
s_prev2=0
coeff=2*cos(2*PI*normalized_frequency);
foreachsample,x
,
s=x
+coeff*s_prev-s_prev2;
s_prev2=s_prev;
s_prev=s;
end
power=s_prev2*s_prev2+s_prev*s_prev-coeff*s_prev2*s_prev;
[
Inordertocomputeasingle
Whenthenumberofdesired
orif,forsomereason,Nisnotanintegralpowerof2whileyousticktoasimpleFFTalgorithmlikethe2-radix
Moreover,theGoertzelalgorithmcanbecomputedassamplescomein,andtheFFTalgorithmmayrequirealargetableofNpre-computedsinesandcosinesinordertobeefficient.
Ifmultiplicationsarenotconsideredasdifficultasadditions,orviceversa,the5/6ratiocanshiftbetweenanythingfrom3/4(additionsdominate)to1/1(multiplicationsdominate).
[
Theterm
Atpresent,DTMFdetectorsaremostoftenimplementedasnumericalalgorithmsoneithergeneralpurposecomputersoronfastdigitalsignalprocessors.Thealgorithmshownbelowisanexampleofsuchadetector.
However,thisalgorithmneedsanadditionalpost-processingsteptocompletelyimplementafunctionalDTMFtonedetector.DTMFtoneburstscanbeasshortas50milli-secondsoraslongasseveralseconds.Thetoneburstcanhavenoiseordropoutswithinitwhichmustbeignored.TheGoertzelalgorithmproducesmultipleoutputs;apost-processingstepneedstosmooththeseoutputsintooneoutputpertoneburst.
Oneadditionalproblemisthatthealgorithmwillsometimesproducespuriousoutputsbecauseofawindowperiodthatisnotcompletelyfilledwithsamples.ImagineaDTMFtoneburstandthenimaginethewindowsuperimposedoverthistoneburst.Obviously,thedetectorisrunningatafixedrateandthetoneburstisnotguaranteedtoarrivealignedwiththetimingofthedetector.Sosomewindowintervalsontheleadingandtrailingedgesofthetoneburstwillnotbeentirelyfilledwithvalidtonesamples.Worse,RC-basedtonegeneratorswilloftenhavevoltagesag/surgerelatedanomaliesattheleadingandtrailingedgesofthetoneburst.Thesealsocancontributetospuriousoutputs.
Itishighlylikelythatthisdetectorwillreportfalseorincorrectresultsattheleadingandtrailingedgesofthetoneburstduetoalackofsufficientvalidsampleswithinthewindow.Inaddition,thetonedetectormustbeabletotoleratetonedropoutswithinthetoneburstandthesecanproduceadditionalfalsereportsduetothesamewindowingeffects.
Thepost-processingsystemcanbeimplementedasastatisticalaggregatorwhichwillexamineaseriesofoutputsofthealgorithmbelow.Thereshouldbeacounterforeachpossibleoutput.Theseallstartoutatzero.Thedetectorstartsproducingoutputsanddependingontheoutput,theappropriatecounterisincremented.Finally,thedetectorstopsgeneratingoutputsforlongenoughthatthetoneburstcanbeconsideredtobeover.ThecounterwiththehighestvaluewinsandshouldbeconsideredtobetheDTMFdigitsignaledbythetoneburst.
WhileitistruethatthereareeightpossiblefrequenciesinaDTMFtone,thealgorithmasoriginallyenteredonthispagewascomputingafewmorefrequenciessoastohelprejectfalsetones(talkoff).Noticethepeaktonecounterloop.Thischeckstoseethatonlytwotonesareactive.Ifmorethanthisarefoundthenthetoneisrejected.
[
#defineSAMPLING_RATE8000
#defineMAX_BINS8
#defineGOERTZEL_N92
intsample_count;
doubleq1[MAX_BINS];
doubleq2[MAX_BINS];
doubler[MAX_BINS];
/*
*coef=2.0*cos((2.0*PI*k)/(float)GOERTZEL_N));
*Wherek=(int)(0.5+((float)GOERTZEL_N*target_freq)/SAMPLING_RATE));
*
*Moresimply:coef=2.0*cos((2.0*PI*target_freq)/SAMPLING_RATE);
*/
doublefreqs[MAX_BINS]=
{
697,
770,
852,
941,
1209,
1336,
1477,
1633
};
doublecoefs[MAX_BINS];
/*----------------------------------------------------------------------------
*calc_coeffs
*----------------------------------------------------------------------------
*Thisiswherewecalculatethecorrectco-efficients.
*/
voidcalc_coeffs()
{
intn;
for(n=0;n<MAX_BINS;n++)
{
coefs[n]=2.0*cos(2.0*3.141592654*freqs[n]/SAMPLING_RATE);
}
}
/*----------------------------------------------------------------------------
*post_testing
*----------------------------------------------------------------------------
*Thisiswherewelookatthebinsanddecideifwehaveavalidsignal.
*/
voidpost_testing()
{
introw,col,see_digit;
intpeak_count,max_index;
doublemaxval,t;
inti;
char*row_col_ascii_codes[4][4]={
{"1","2","3","A"},
{"4","5","6","B"},
{"7","8","9","C"},
{"*","0","#","D"}};
/*Findthelargestintherowgroup.*/
row=0;
maxval=0.0;
for(i=0;i<4;i++)
{
if(r>maxval)
{
maxval=r[i];
row=i;
}
}
[i]/*Findthelargestinthecolumngroup.*/
col=4;
maxval=0.0;
for(i=4;i<8;i++)
{
if(r>maxval)
{
maxval=r[i];
col=i;
}
}
[i]/*Checkforminimumenergy*/
if(r[row]<4.0e5)/*2.0e5...1.0e8nochange*/
{
/*energynothighenough*/
}
elseif(r[col]<4.0e5)
{
/*energynothighenough*/
}
else
{
see_digit=TRUE;
/*Twistcheck
*CEPT=>twist<6dB
*AT&T=>forwardtwist<4dBandreversetwist<8dB
*-ndB<10log10(v1/v2),wherev1<v2
*-4dB<10log10(v1/v2)
*-0.4<log10(v1/v2)
*0.398<v1/v2
*0.398*v2<v1
*/
if(r[col]>r[row])
{
/*Normaltwist*/
max_index=col;
if(r[row]<(r[col]*0.398))/*twist>4dB,error*/
see_digit=FALSE;
}
else/*if(r[row]>r[col])*/
{
/*Reversetwist*/
max_index=row;
if(r[col]<(r[row]*0.158))/*twist>8db,error*/
see_digit=FALSE;
}
/*Signaltonoisetest
*AT&Tstatesthatthenoisemustbe16dBdownfromthesignal.
*Herewecountthenumberofsignalsabovethethresholdand
*thereoughttobeonlytwo.
*/
if(r[max_index]>1.0e9)
t=r[max_index]*0.158;
else
t=r[max_index]*0.010;
peak_count=0;
for(i=0;i<8;i++)
{
if(r>t)
peak_count++;
}
if(peak_count>2)
see_digit=FALSE;
if(see_digit)
{
printf("%s",row_col_ascii_codes[row][col-4]);
fflush(stdout);
}
}
}
/*----------------------------------------------------------------------------[i]
*goertzel
*----------------------------------------------------------------------------
*/
voidgoertzel(intsample)
{
doubleq0;
ui32i;
sample_count++;
/*q1[0]=q2[0]=0;*/
for(i=0;i<MAX_BINS;i++)
{
q0=coefs*q1[i]-q2[i]+sample;
q2[i]=q1[i];
q1[i]=q0;
}
if(sample_count==GOERTZEL_N)
{
for(i=0;i<MAX_BINS;i++)
{
r[i]=(q1[i]*q1[i])+(q2[i]*q2[i])-(coefs[i]*q1[i]*q2[i]);
q1[i]=0.0;
q2[i]=0.0;
}
post_testing();
sample_count=0;
}
}
[
[
TheGoertzelAlgorithm
TheGoertzelalgorithmcanperformtonedetectionusingmuchlessCPUhorsepowerthantheFastFourierTransform,butmanyengineershaveneverheardofit.Thisarticleattemptstochangethat.
MostengineersarefamiliarwiththeFastFourierTransform(FFT)andwouldhavelittletroubleusinga"canned"FFTroutinetodetectoneormoretonesinanaudiosignal.Whatmanydon'tknow,however,isthatifyouonlyneedtodetectafewfrequencies,amuchfastermethodisavailable.It'scalledtheGoertzelalgorithm.
Tonedetection
Manyapplicationsrequiretonedetection,suchas:
DTMF(touchtone)decoding
Callprogress(dialtone,busy,andsoon)decoding
Frequencyresponsemeasurements(sendingatonewhilesimultaneouslyreadingbacktheresult)-ifyoudothisforarangeoffrequencies,theresultingfrequencyresponsecurvecanbeinformative.Forexample,thefrequencyresponsecurveofatelephonelinetellsyouifanyloadcoils(inductors)arepresentonthatline.
AlthoughdedicatedICsexistfortheapplicationsabove,implementingthesefunctionsinsoftwarecostsless.Unfortunately,manyembeddedsystemsdon'thavethehorsepowertoperformcontinuousreal-timeFFTs.That'swheretheGoertzelalgorithmcomesin.
Inthisarticle,IdescribewhatIcallabasicGoertzelandanoptimizedGoertzel.
ThebasicGoertzelgivesyourealandimaginaryfrequencycomponentsasaregularDiscreteFourierTransform(DFT)orFFTwould.Ifyouneedthem,magnitudeandphasecanthenbecomputedfromthereal/imaginarypair.
TheoptimizedGoertzelisevenfaster(andsimpler)thanthebasicGoertzel,butdoesn'tgiveyoutherealandimaginaryfrequencycomponents.Instead,itgivesyoutherelativemagnitudesquared.Youcantakethesquarerootofthisresulttogettherelativemagnitude(ifneeded),butthere'snowaytoobtainthephase.
Inthisshortarticle,Iwon'ttrytoexplainthetheoreticalbackgroundofthealgorithm.Idogivesomelinksattheendwhereyoucanfindmoredetailedexplanations.Icantellyouthatthealgorithmworkswell,havinguseditinallofthetonedetectionapplicationspreviouslylisted(andothers).
AbasicGoertzel
Firstaquickoverviewofthealgorithm:someintermediateprocessingisdoneineverysample.TheactualtonedetectionoccurseveryNthsample.(I'lltalkmoreaboutNinaminute.)
AswiththeFFT,youworkwithblocksofsamples.However,thatdoesn'tmeanyouhavetoprocessthedatainblocks.Thenumericalprocessingisshortenoughtobedoneintheveryinterruptserviceroutine(ISR)thatisgatheringthesamples(ifyou'regettinganinterruptpersample).Or,ifyou'regettingbuffersofsamples,youcangoaheadandprocessthemabatchatatime.
BeforeyoucandotheactualGoertzel,youmustdosomepreliminarycalculations:
Decideonthesamplingrate.
Choosetheblocksize,N.
Precomputeonecosineandonesineterm.
Precomputeonecoefficient.
Thesecanallbeprecomputedonceandthenhardcodedinyourprogram,savingRAMandROMspace;oryoucancomputethemon-the-fly.
Samplingrate
Yoursamplingratemayalreadybedeterminedbytheapplication.Forexample,intelecomapplications,it'scommontouseasamplingrateof8kHz(8,000samplespersecond).Alternatively,youranalog-to-digitalconverter(orCODEC)mayberunningfromanexternalclockorcrystaloverwhichyouhavenocontrol.
Ifyoucanchoosethesamplingrate,theusualNyquistrulesapply:thesamplingratewillhavetobeatleasttwiceyourhighestfrequencyofinterest.Isay"atleast"becauseifyouaredetectingmultiplefrequencies,it'spossiblethatanevenhighersamplingfrequencywillgivebetterresults.Whatyoureallywantisforeveryfrequencyofinteresttobeanintegerfactorofthesamplingrate.
Blocksize
GoertzelblocksizeNislikethenumberofpointsinanequivalentFFT.Itcontrolsthefrequencyresolution(alsocalledbinwidth).Forexample,ifyoursamplingrateis8kHzandNis100samples,thenyourbinwidthis80Hz.
ThiswouldsteeryoutowardsmakingNashighaspossible,togetthehighestfrequencyresolution.ThecatchisthatthehigherNgets,thelongerittakestodetecteachtone,simplybecauseyouhavetowaitlongerforallthesamplestocomein.Forexample,at8kHzsampling,itwilltake100msfor800samplestobeaccumulated.Ifyou'retryingtodetecttonesofshortduration,youwillhavetousecompatiblevaluesofN.
ThethirdfactorinfluencingyourchoiceofNistherelationshipbetweenthesamplingrateandthetargetfrequencies.Ideallyyouwantthefrequenciestobecenteredintheirrespectivebins.Inotherwords,youwantthetargetfrequenciestobeintegermultiplesofsample_rate/N.
Thegoodnewsisthat,unliketheFFT,Ndoesn'thavetobeapoweroftwo.
Precomputedconstants
Onceyou'veselectedyoursamplingrateandblocksize,it'sasimplefive-stepprocesstocomputetheconstantsyou'llneedduringprocessing:
[i]w
cosine=cosw
sine=sinw
coeff=2*cosine
Fortheper-sampleprocessingyou'regoingtoneedthreevariables.Let'scallthemQ0,Q1,andQ2.
Q1isjustthevalueofQ0lasttime.Q2isjustthevalueofQ0twotimesago(orQ1lasttime).
Q1andQ2mustbeinitializedtozeroatthebeginningofeachblockofsamples.Foreverysample,youneedtorunthefollowingthreeequations:
Q0=coeff*Q1-Q2+sample
Q2=Q1
Q1=Q0
Afterrunningtheper-sampleequationsNtimes,it'stimetoseeifthetoneispresentornot.
real=(Q1-Q2*cosine)
imag=(Q2*sine)
magnitude2=real2+imag2
Asimplethresholdtestofthemagnitudewilltellyouifthetonewaspresentornot.ResetQ2andQ1tozeroandstartthenextblock.
AnoptimizedGoertzel
TheoptimizedGoertzelrequireslesscomputationthanthebasicone,attheexpenseofphaseinformation.
Theper-sampleprocessingisthesame,buttheendofblockprocessingisdifferent.Insteadofcomputingrealandimaginarycomponents,andthenconvertingthoseintorelativemagnitudesquared,youdirectlycomputethefollowing:
magnitude2=Q12+Q22-Q1*Q2*coeff
ThisistheformofGoertzelI'veusedmostoften,anditwasthefirstoneIlearnedabout.
Pullingitalltogether
Theprogramdoestwodemonstrations.Inthefirstone,bothformsoftheGoertzelalgorithmareusedtocomputerelativemagnitudesquaredandrelativemagnitudeforthreedifferentsynthesizedsignals:onebelowthetarget_frequency,oneequaltothetarget_frequency,andoneabovethetarget_frequency.
You'llnoticethattheresultsarenearlyidentical,regardlessofwhichformoftheGoertzelalgorithmisused.You'llalsonoticethatthedetectorvaluespeaknearthetargetfrequency.
Intheseconddemonstration,asimulatedfrequencysweepisrun,andtheresultsofjustthebasicGoertzelareshown.Again,you'llnoticeaclearpeakinthedetectoroutputnearthetargetfrequency.Figure1showstheoutputofthecodein
Figure1:Demooutput
ForSAMPLING_RATE=8000.000000N=205andFREQUENCY=941.000000,k=24andcoeff=1.482867
Fortestfrequency691.000000:
real=-360.392059imag=-45.871609
Relativemagnitudesquared=131986.640625
Relativemagnitude=363.299652
Relativemagnitudesquared=131986.640625
Relativemagnitude=363.299652
Fortestfrequency941.000000:
real=-3727.528076imag=-9286.238281
Relativemagnitudesquared=100128688.000000
Relativemagnitude=10006.432617
Relativemagnitudesquared=100128680.000000
Relativemagnitude=10006.431641
Fortestfrequency1191.000000:
real=424.038116imag=-346.308716
Relativemagnitudesquared=299738.062500
Relativemagnitude=547.483398
Relativemagnitudesquared=299738.062500
Relativemagnitude=547.483398
| Freq= | 641.0 | relmag^2= | 146697.87500 | relmag= | 383.01160 |
| Freq= | 656.0 | relmag^2= | 63684.62109 | relmag= | 252.35812 |
| Freq= | 671.0 | relmag^2= | 96753.92188 | relmag= | 311.05292 |
| Freq= | 686.0 | relmag^2= | 166669.90625 | relmag= | 408.25226 |
| Freq= | 701.0 | relmag^2= | 5414.02002 | relmag= | 73.58002 |
| Freq= | 716.0 | relmag^2= | 258318.37500 | relmag= | 508.25031 |
| Freq= | 731.0 | relmag^2= | 178329.68750 | relmag= | 422.29099 |
| Freq= | 746.0 | relmag^2= | 71271.88281 | relmag= | 266.96796 |
| Freq= | 761.0 | relmag^2= | 437814.90625 | relmag= | 661.67584 |
| Freq= | 776.0 | relmag^2= | 81901.81250 | relmag= | 286.18494 |
| Freq= | 791.0 | relmag^2= | 468060.31250 | relmag= | 684.14935 |
| Freq= | 806.0 | relmag^2= | 623345.56250 | relmag= | 789.52234 |
| Freq= | 821.0 | relmag^2= | 18701.58984 | relmag= | 136.75375 |
| Freq= | 836.0 | relmag^2= | 1434181.62500 | relmag= | 1197.57324 |
| Freq= | 851.0 | relmag^2= | 694211.75000 | relmag= | 833.19373 |
| Freq= | 866.0 | relmag^2= | 1120359.50000 | relmag= | 1058.47034 |
| Freq= | 881.0 | relmag^2= | 4626623.00000 | relmag= | 2150.95874 |
| Freq= | 896.0 | relmag^2= | 160420.43750 | relmag= | 400.52521 |
| Freq= | 911.0 | relmag^2= | 19374364.00000 | relmag= | 4401.63184 |
| Freq= | 926.0 | relmag^2= | 81229848.00000 | relmag= | 9012.76074 |
| Freq= | 941.0 | relmag^2= | 100128688.00000 | relmag= | 10006.43262 |
| Freq= | 956.0 | relmag^2= | 43694608.00000 | relmag= | 6610.18994 |
| Freq= | 971.0 | relmag^2= | 1793435.75000 | relmag= | 1339.19226 |
| Freq= | 986.0 | relmag^2= | 3519388.50000 | relmag= | 1876.00330 |
| Freq= | 1001.0 | relmag^2= | 3318844.50000 | relmag= | 1821.76965 |
| Freq= | 1016.0 | relmag^2= | 27707.98828 | relmag= | 166.45717 |
| Freq= | 1031.0 | relmag^2= | 1749922.62500 | relmag= | 1322.84644 |
| Freq= | 1046.0 | relmag^2= | 478859.28125 | relmag= | 691.99658 |
| Freq= | 1061.0 | relmag^2= | 284255.81250 | relmag= | 533.15643 |
| Freq= | 1076.0 | relmag^2= | 898392.93750 | relmag= | 947.83594 |
| Freq= | 1091.0 | relmag^2= | 11303.36035 | relmag= | 106.31726 |
| Freq= | 1106.0 | relmag^2= | 420975.65625 | relmag= | 648.82635 |
| Freq= | 1121.0 | relmag^2= | 325753.78125 | relmag= | 570.74841 |
| Freq= | 1136.0 | relmag^2= | 36595.78906 | relmag= | 191.30026 |
| Freq= | 1151.0 | relmag^2= | 410926.06250 | relmag= | 641.03516 |
| Freq= | 1166.0 | relmag^2= | 45246.58594 | relmag= | 212.71245 |
| Freq= | 1181.0 | relmag^2= | 119967.59375 | relmag= | 346.36337 |
| Freq= | 1196.0 | relmag^2= | 250361.39062 | relmag= | 500.36127 |
| Freq= | 1211.0 | relmag^2= | 1758.44263 | relmag= | 41.93379 |
| Freq= | 1226.0 | relmag^2= | 190195.57812 | relmag= | 436.11417 |
| Freq= | 1241.0 | relmag^2= | 74192.23438 | relmag= | 272.38251 |
Asyoucansee,theGoertzelalgorithmdeservestobeaddedtoyoursignalprocessingtoolbox.
KevinBankshasbeendevelopingembeddedsystemsfor19years,asaconsultantandasanemployeeatcompaniesincludingSCI,TxPORT,DiscoveryCom,andNokia.Currentlyheisbackinconsultingmode.Hecanbereachedat
References
Herearesomelinkstootherresourcesyoumayfinduseful:
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 