Goertzel algoritm

Goertzel -algoritmen är en teknik inom digital signalbehandling (DSP) för effektiv utvärdering av de individuella termerna för den diskreta Fouriertransformen (DFT). Det är användbart i vissa praktiska tillämpningar, såsom igenkänning av dual-tone multi-frequency signaling (DTMF) toner som produceras av tryckknapparna på knappsatsen på en traditionell analog telefon . Algoritmen beskrevs första gången av Gerald Goertzel 1958.

Precis som DFT analyserar Goertzel-algoritmen en valbar frekvenskomponent från en diskret signal . Till skillnad från direkta DFT-beräkningar tillämpar Goertzel-algoritmen en enda reellt värderad koefficient vid varje iteration, med hjälp av realvärderad aritmetik för realvärderade indatasekvenser. För att täcka ett fullt spektrum (förutom när den används för kontinuerlig dataström där koefficienter återanvänds för efterföljande beräkningar, som har beräkningskomplexitet som motsvarar glidande DFT ), har Goertzel-algoritmen en högre komplexitetsordning än snabba Fourier-transformeringsalgoritmer (FFT), men för att beräkna ett litet antal utvalda frekvenskomponenter är det mer numeriskt effektivt. Den enkla strukturen hos Goertzel-algoritmen gör den väl lämpad för små processorer och inbyggda applikationer.

Goertzel-algoritmen kan också användas "omvänt" som en sinusformad syntesfunktion, som endast kräver 1 multiplikation och 1 subtraktion per genererat prov.

Algoritmen

Huvudberäkningen i Goertzel-algoritmen har formen av ett digitalt filter , och av denna anledning kallas algoritmen ofta för ett Goertzel-filter . Filtret arbetar på en ingångssekvens $x[n]$ i en kaskad av två steg med en parameter $\omega _{0}$ , vilket ger frekvensen som ska analyseras, normaliserad till radianer per prov.

Det första steget beräknar en mellansekvens, $s[n]$ :

$s[n]=x[n]+2\cos(\omega _{0})s[n-1]-s[n-2].$

 

 

 

 

(1)

Det andra steget tillämpar följande filter på $s[n]$ , och producerar utdatasekvensen $y[n]$ :

$y[n]=s[n]-e^{-j\omega _{0}}s[n-1].$

 

 

 

 

(2)

Det första filtersteget kan observeras vara ett andra ordningens IIR-filter med en direktformad struktur. Denna speciella struktur har egenskapen att dess interna tillståndsvariabler är lika med tidigare utdatavärden från det steget. Indatavärden $x[n]$ för $n<0$ antas alla vara lika med 0. För att fastställa det initiala filtertillståndet så att utvärderingen kan börja vid exempel $x[0]$ , filtertillstånden tilldelas initiala värden $s[-2]=s[-1]=0$ . För att undvika aliasrisker är frekvensen $\omega _{0}$ ofta begränsad till intervallet 0 till π (se Nyquist–Shannons samplingssats ) ; att använda ett värde utanför detta intervall är inte meningslöst, men det motsvarar att använda en aliaserad frekvens inom detta intervall, eftersom exponentialfunktionen är periodisk med en period på 2π i ω {\ $\omega _{0}}$ .

Andrastegsfiltret kan observeras vara ett FIR-filter , eftersom dess beräkningar inte använder någon av dess tidigare utdata.

Z- transformmetoder kan användas för att studera egenskaperna hos filterkaskaden. Z-transformen för det första filtersteget som ges i ekvation (1) är

${\begin{aligned}{\frac {S(z)}{X(z)}}&={\frac {1}{1-2\cos(\omega _{0})z^{ -1}+z^{-2}}}\\&={\frac {1}{(1-e^{+j\omega _{0}}z^{-1})(1-e^ {-j\omega _{0}}z^{-1})}}.\end{aligned}}$

 

 

 

 

(3)

Z-transformen för det andra filtersteget som ges i ekvation (2) är

${\frac {Y(z)}{S(z)}}=1-e^{-j\omega _{0}}z^{-1}.$

 

 

 

 

(4)

Den kombinerade överföringsfunktionen för kaskaden av de två filterstegen är då

${\begin{aligned}{\frac {S(z)}{X(z)}}{\frac {Y(z)}{S(z)}}={\frac {Y(z) }{X(z)}}&={\frac {(1-e^{-j\omega _{0}}z^{-1})}{(1-e^{+j\omega _{ 0}}z^{-1})(1-e^{-j\omega _{0}}z^{-1})}}\\&={\frac {1}{1-e^{ +j\omega _{0}}z^{-1}}}.\end{aligned}}$

 

 

 

 

(5)

Detta kan transformeras tillbaka till en ekvivalent tidsdomänsekvens, och termerna rullas tillbaka till den första inmatningstermen vid index $n=0$ :

${\begin{aligned}y[n]&=x[n]+e^{+j\omega _{0 }}y[n-1]\\&=\summa _{k=-\infty }^{n}x[k]e^{+j\omega _{0}(nk)}\\&=e ^{j\omega _{0}n}\summa _{k=0}^{n}x[k]e^{-j\omega _{0}k}\qquad {\text{sedan }}\ forall k<0,x[k]=0.\end{aligned}}$

 

 

 

 

(6)

Numerisk stabilitet

Det kan observeras att polerna för filtrets Z-transform är belägna vid $e^{+j\omega _{0}}$ och $e^{-j\omega _{0}}$ , på en cirkel med enhetsradie centrerad på ursprunget för det komplexa Z-transformplanet. Denna egenskap indikerar att filterprocessen är marginellt stabil och sårbar för ackumulering av numeriska fel när den beräknas med aritmetik med låg precision och långa ingångssekvenser. En numeriskt stabil version föreslogs av Christian Reinsch.

DFT-beräkningar

För det viktiga fallet med att beräkna en DFT-term tillämpas följande speciella begränsningar.

• Filtreringen slutar vid index $n=N$ , där $N$ är antalet termer i inmatningssekvensen för DFT.
• De frekvenser som valts för Goertzel-analysen är begränsade till den speciella formen

$\omega _{0}=2\pi {\frac {k}{N}}.$

 

 

 

 

(7)

• Indexnumret $k$ som indikerar "frekvensfacket" för DFT väljs från uppsättningen av indexnummer

$k\in \{0,1,2,...,N-1\}.$

 

 

 

 

(8)

Genom att göra dessa substitutioner till ekvation (6) och observera att termen $e^{+j2\pi k}=1$ , ekvation (6) tar sedan följande form:

$y[N]=\sum _{n=0}^{N}x[n]e^{-j2\pi {\frac {nk}{N}}}.$

 

 

 

 

(9)

Vi kan observera att den högra sidan av ekvation (9) är extremt lik den definierande formeln för DFT-term ${\displaystyle X[k]} ,$ DFT-termen för indexnummer $k$ , men inte exakt samma. Summan som visas i ekvation (9) kräver $N+1$ inmatningstermer, men endast $N$ indatatermer är tillgängliga vid utvärdering av en DFT. En enkel men oelegant lösning är att utöka inmatningssekvensen $x[n]$ med ytterligare ett artificiellt värde $x[N]=0$ . Vi kan se från ekvation (9) att den matematiska effekten på slutresultatet är densamma som att ta bort termen $x[N]$ från summeringen, och därmed leverera det avsedda DFT-värdet.

Det finns dock ett mer elegant tillvägagångssätt som undviker det extra filterpasset. Från ekvation (1) kan vi notera att när den utökade inmatningstermen $x[N]=0$ används i det sista steget,

$s[N]=2\cos(\omega _{0})s[N-1]-s[N-2].$

 

 

 

 

(10)

Således kan algoritmen slutföras enligt följande:

• avsluta IIR-filtret efter bearbetning av inmatningsterm $x[N-1]$ ,
• tillämpa ekvation (10) för att konstruera $s[N]$ från de tidigare utgångarna $s[N-1]$ och $s [N-2]$ ,
• tillämpa ekvation (2) med det beräknade ${\displaystyle s[N]}-$ värdet och med $s[N-1]$ producerat av den slutliga direkta beräkningen av filtret.

De två sista matematiska operationerna förenklas genom att kombinera dem algebraiskt:

${\begin{aligned}y[N]&=s[N]-e^{-j2\pi {\frac {k}{N}}}s[N-1]\\&=(2 \cos(\omega _{0})s[N-1]-s[N-2])-e^{-j2\pi {\frac {k}{N}}}s[N-1]\ \&=e^{j2\pi {\frac {k}{N}}}s[N-1]-s[N-2].\end{aligned}}$

 

 

 

 

(11)

Observera att om du stoppar filteruppdateringarna vid term $N-1$ och omedelbart tillämpar ekvation (2) snarare än ekvation (11), missas de slutliga filtertillståndsuppdateringarna, vilket ger ett resultat med felaktig fas.

Den speciella filtreringsstrukturen som valts för Goertzel-algoritmen är nyckeln till dess effektiva DFT-beräkningar. Vi kan observera att endast ett utdatavärde $y[N]$ används för att beräkna DFT, så beräkningar för alla andra utmatningstermer utelämnas. Eftersom FIR-filtret inte beräknas, kan IIR-stegsberäkningarna $s[0],s[1]$ etc. kasseras omedelbart efter uppdatering av det första stegets interna tillstånd.

Detta tycks lämna en paradox: för att slutföra algoritmen måste FIR-filtersteget utvärderas en gång med användning av de två sista utgångarna från IIR-filtersteget, medan för beräkningseffektivitet kasserar IIR-filtrets iteration sina utvärden. Det är här egenskaperna för filterstrukturen i direkt form tillämpas. De två interna tillståndsvariablerna för IIR-filtret tillhandahåller de två sista värdena på IIR-filtrets utgång, vilka är termerna som krävs för att utvärdera FIR-filtersteget.

Ansökningar

Effektspektrumtermer

Genom att undersöka ekvation (6), ett slutligt IIR-filterpass för att beräkna termen $y[N]$ med hjälp av ett tilläggsvärde $x[N]=0$ tillämpar en komplex multiplikator av magnitud 1 till föregående term $y[N-1]$ . Följaktligen $y[N]$ och $y[N-1]$ ekvivalent signaleffekt. Det är lika giltigt att tillämpa ekvation (11) och beräkna signaleffekten från term $y[N]$ eller att tillämpa ekvation (2) och beräkna signaleffekten från term ${ \displaystyle y[N-1]}$ . Båda fallen leder till följande uttryck för signaleffekten som representeras av DFT-term $X[k]$ :

${\begin{aligned}X[k]\,X'[k]&=y[N]\,y'[N]=y[N-1]\,y'[N-1]\ \&=s^{2}[N-1]+s^{2}[N-2]-2\cos \left(2\pi {\frac {k}{N}}\höger)\,s [N-1]\,s[N-2].\end{aligned}}$

 

 

 

 

(12)

I pseudokoden nedan lagras indata med verkligt värde i arrayen x och variablerna sprev och sprev2 lagrar temporärt utdatahistorik från IIR-filtret. Nterms är antalet sampel i arrayen, och Kterm motsvarar frekvensen av intresse, multiplicerat med samplingsperioden.

 Ntermer definierade här Kterm vald här ω = 2 × π × Kterm / Nterms; koeff := 2 × cos(ω) sprev := 0 sprev2 := 0   för varje  index  n  i intervallet 0 till Nterms-1  do  s := x[n] + coeff × sprev - sprev2 sprev2 := sprev sprev := s  sluteffekt  := sprev  2  + sprev2  2  - (koeff × sprev × sprev2) 

Det är möjligt att organisera beräkningarna så att inkommande prover levereras var för sig till ett mjukvaruobjekt som upprätthåller filtertillståndet mellan uppdateringar, med det slutliga effektresultatet åtkomligt efter att den andra behandlingen är klar.

Enskild DFT-term med realvärderad aritmetik

Fallet med realvärderade indata uppstår ofta, särskilt i inbyggda system där ingångsströmmarna härrör från direkta mätningar av fysiska processer. När indata är realvärderade kan filterinterna tillståndsvariabler sprev och sprev2 också observeras vara verkligt värderade, följaktligen krävs ingen komplex aritmetik i det första IIR-steget. Att optimera för realvärderad aritmetik är vanligtvis så enkelt som att tillämpa lämpliga realvärderade datatyper för variablerna.

Efter att beräkningarna med inmatningsterm $x[N-1]$ och filteriterationer har avslutats, måste ekvation (11) tillämpas för att utvärdera DFT-termen. Den slutliga beräkningen använder komplext värderad aritmetik, men denna kan omvandlas till realvärderad aritmetik genom att separera reella och imaginära termer:

${\begin{aligned}c_{r}&=\cos(2\pi {\tfrac {k}{N}}),\\c_{i}&=\sin(2\pi {\tfrac {k}{N}}),\\y[N]&=c_{r}s[N-1]-s[N-2]+jc_{i}s[N-1].\end{aligned }}$

 

 

 

 

(13)

Jämfört med effektspektrumapplikationen är den enda skillnaden beräkningen som används för att avsluta:

(Samma IIR-filterberäkningar som i signaleffektimplementeringen) XKreal = sprev * cr - sprev2; XKimag = sprev * ci;

Fasdetektering

Den här applikationen kräver samma utvärdering av DFT-term ${\displaystyle X[k]} ,$ som diskuterades i föregående avsnitt, med användning av en indataström med verkligt värde eller komplext värde. Sedan kan signalfasen utvärderas som

$\phi =\tan ^{-1}{\frac {\Im (X[k])}{\Re (X[k])}},$

 

 

 

 

(14)

vidta lämpliga försiktighetsåtgärder för singulariteter, kvadranter och så vidare vid beräkning av den inversa tangentfunktionen.

Komplexa signaler i verklig aritmetik

Eftersom komplexa signaler bryts ner linjärt till verkliga och imaginära delar, kan Goertzel-algoritmen beräknas i verklig aritmetik separat över sekvensen av reella delar, vilket ger $y_{\text{r}}[n]$ , och över sekvensen av imaginära delar, vilket ger $y_{\text{i}}[n]$ . Därefter kan de två komplext värderade delresultaten kombineras igen:

$y[n]=y_{\text{r}}[n]+jy_{\text{i}}[n].$

 

 

 

 

(15)

Beräkningskomplexitet

• Enligt beräkningskomplexitetsteori , beräkna en uppsättning $M$ DFT-termer med $M$ -applikationer av Goertzel-algoritmen på en datamängd med $N$ -värden med en "kostnad per operation" på $K$ har komplexitet $O(KNM)$ .

För att beräkna en enda DFT bin $X(f)$ för en komplex ingångssekvens med längden $N$ kräver Goertzel-algoritmen $2N$ multiplikationer och $4\ N$ additioner/subtraktioner inom slingan, samt 4 multiplikationer och 4 slutliga additioner/subtraktioner, för totalt $2N+4$ multiplikationer och $4N+4$ tillägg/subtraktioner. Detta upprepas för var och en av $M$ -frekvenserna.

• Om man däremot använder en FFT på en datamängd med $N$ -värden har komplexiteten $O(KN\log N)$ .

Detta är svårare att tillämpa direkt eftersom det beror på vilken FFT-algoritm som används, men ett typiskt exempel är en radix-2 FFT, som kräver $2\log _{2}(N)$ multiplikationer och $3\log _{2}(N)$ additioner/subtraktioner per DFT- fack, för var och en av de $N$ -fackarna.

I komplexitetsordningsuttrycken, när antalet beräknade termer $M$ är mindre än $\log N$ , är fördelen med Goertzel-algoritmen tydlig. Men eftersom FFT-koden är jämförelsevis komplex är faktorn "kostnad per arbetsenhet" $K$ ofta större för en FFT, och den praktiska fördelen gynnar Goertzel-algoritmen även för $M$ flera gånger större än $\log _{2}(N)$ .

Som en tumregel för att avgöra om en radix-2 FFT eller en Goertzel-algoritm är mer effektiv, justera antalet termer $N$ i datauppsättningen uppåt till närmaste exakta potens av 2, och kalla detta $N_{2}$ och Goertzel-algoritmen är sannolikt snabbare om

$M\leq {\frac {5N_{2}}{6N}}\log _{2}(N_{2})$

FFT-implementeringar och bearbetningsplattformar har en betydande inverkan på den relativa prestandan. Vissa FFT-implementeringar utför interna komplexa talberäkningar för att generera koefficienter i farten, vilket avsevärt ökar deras "kostnad K per arbetsenhet." FFT- och DFT-algoritmer kan använda tabeller med förberäknade koefficientvärden för bättre numerisk effektivitet, men detta kräver mer åtkomst till koefficientvärden buffrade i externt minne, vilket kan leda till ökad cachekonflikt som motverkar en del av de numeriska fördelarna.

Båda algoritmerna får ungefär en faktor 2 effektivitet när man använder realvärderade snarare än komplexa indata. Dessa vinster är dock naturliga för Goertzel-algoritmen men kommer inte att uppnås för FFT utan att använda vissa algoritmvarianter [ specialiserat på att transformera verkligt värderade data .

Se även

• Bluesteins FFT-algoritm (chirp-Z)
• Frequency Shift Keying (FSK)
• Fasskiftningsnyckel (PSK)

Vidare läsning

• Proakis, JG; Manolakis, DG (1996), Digital Signal Processing: Principles, Algorithms, and Applications , Upper Saddle River, NJ: Prentice Hall, s. 480–481, Bibcode : 1996dspp.book.....P

externa länkar

• Goertzel Algorithm at the Wayback Machine (arkiverad 2018-06-28)
• En DSP-algoritm för frekvensanalys
• Goertzel-algoritmen av Kevin Banks