Goertzel-Algorithmus und Signed Fractional Format

[Tony P. (micropower)]

Hallo!

Um (breitbandige und kontinuierliche) akustischen Signalen mit einem 
Mikrophon Array zu lokalisieren filtere ich zurzeit eine ausgewählte 
Frequenz (Wellenlänge = 2 * Mikrophonabstand) und berechne anschließend 
den cross-correlation-coefficient (ccc) zwischen den Signalen.

Diese Methode ist sehr rechenintensiv, weiterhin unterscheidet das 
CCC-Ergebnis nicht zwischen Rechts und Links. Eine Unterscheidung 
zwischen Rechts und Links ist nur möglich, wenn ich die Signale 
verschiebe und den CCC erneut berechne.

Nun möchte ich mit Goertzel-Algorithmus mein Glück versuchen und die 
Phasenverschiebung berechnen.

Meine Versuche sind aber fehlgeschlagen. Vielleicht habe ich die Formeln 
aus http://en.wikipedia.org/wiki/Goertzel_algorithm falsch 
interpretiert.

Zurzeit habe ich meine Samples (Input: x(n)) in Signed Fractional Format 
und 8 Samples pro Welle. Die Frequenz, die mich interessiert ist 1429 Hz 
und die Samplefrequenz ist 11432 Hz

Output = s(n) = x(n) + Coeff * s(n − 1) − s(n − 2)
Coeff = 2cos(2πω) = 2 * 0,707106781
Π = Kreiszahl = 3,141592654
ω = Cycles per Sample = 0,125

Ist meine Berechnung überhaupt richtig?
Der dsPIC Akkumulator ist nach einem Block von 64 Samples überflutet, 
was mache ich falsch?

Danke Euch!

[Tony P. (micropower)]

Sorry, der Akkumulator wurde wegen eines Programmierfehlers überflutet.

Bitte, nicht erschlagen, ich bin ein Anfänger zu DSPs und das ist ein 
privates Projekt (Lernprojekt).

Also, die Formel funktioniert, …
Trotzdem habe ich damit aufgegeben.

"imag" und "real" sind damit schnell berechet, aber, um die Phase zu 
bekommen, muss ich ja jedes Mal arctan(imag/real) berechnen. Mit 
Nährungsformeln ist das natürlich machbar, das ist aber auch 
Rechenintensiv. Diese Methode ist wirklich effektiv, um die Magnitude^2 
für eine Frequenz schnell zu berechnen, aber nicht die Phasen.

Vielen Dank

[Detlef _. (detlef_a)]

Den arctan berechnest Du einfach und schnell mit einem 
Cordic-Algorithmus. Der dreht den komplexen Zeiger und summiert die 
Drehwinkel auf.

Für GCC/WINAVR Code anbei.

Math rulez!
Cheers
Detlef


/************************************************************/

static int32_t wxx[16]={0x40000000,0x25c80a3b,0x13f670b7, 0xa2223a8,
                   0x5161a86, 0x28bafc3, 0x145ec3d,  0xa2f8aa,
                    0x517ca7,  0x28be5d,  0x145f30,   0xa2f98,
                     0x517cc,   0x28be6,   0x145f3,    0xa2fa};

/**************************************************************/
uint16_t IntAtan2(int16_t im, int16_t re)
/**************************************************************/
{ uint32_t uire,uiim;
  uint32_t ww,bla,q,k;


q=0;

//DebugPrintf(" yyy%d \t %d \n",(int32_t)im,(int32_t)re);

if((re==-32768)||(im==-32768))
    {re>>=1;im>>=1;}
//Dreifaltigkeit !
if(im<0 ){im=-im; re=-re;     q+=4;}
if(re<0 ){k =-re; re=im; im=k;q+=2;}
if(re<im){k = re; re=im; im=k;q+=1;}

uire=(uint32_t)re;
uiim=(uint32_t)im;
ww=0;


    if(!(uire|uiim)) return ww;

    while(!((uire|uiim)&(int32_t)0x40000000)){uire<<=1;uiim<<=1;}

    for(k=0;k<16;k++){
        bla=(uiim<<k);
        if((uint32_t)bla<uire){   // zu klein
           uire>>=1;uiim>>=1;
        }else{
           bla -= uire;
           uire=((uire<<k)+uiim)>>(k+1); //Realteil drehen
           uiim=bla>>(k+1);
           ww+=wxx[k];
        }
        if(!uiim) break;
    }

k=ww>>1;
if(q&1) k=((int32_t)1<<30)-k;
k=(((k>>2)+((q&6)<<27))>>14);

return k ;

}

[Tony P. (micropower)]

Vielen Dank für den Hinweis und für den Code. Dieser Cordic-Algorithmus 
ist eine elegante Lösung. Die einzige Lösung, die mir eingefallen war, 
war eine Tabelle mit einigen vorkalkulierten Ergebnissen. Meine Lösung 
habe ich aber nicht verwirklicht gehabt, weil ich mit den komplexen 
Zahlen aufgegeben habe.

Irgendwie glaube ich nicht mehr, dass diese Aufgabe mit 
Goertzel-Algorithmus (oder auch mit FFT und Co.) effektiver gelöst 
werden kann. Trotzdem denke ich immer wieder drüber nach ;-)

Die Phase, die aus "imag" und "real" berechnet wird ist zwar von –Π bis 
+Π, diese Phase hängt aber davon ab, wie die Signale im Sampleblock 
zentriert sind. Sind die Signale unglücklich zentriert, dann bekommt man 
einfach ein falsches Ergebnis.

Vielleicht gibt es für alle diese Probleme eine Lösung? Ich kenne aber 
kein Beispielprojekt.

[Detlef _. (detlef_a)]

>>diese Phase hängt aber davon ab, wie die Signale im Sampleblock zentriert 
>>sind.

So ist es. Aber Du hast ja mehrere Mikrophone, also berechnest Du die 
Phasendifferenz von zweien. Dein erster versuch mit der Korrelation war 
doch auch nicht schlecht. Das bringt Dir zwar keine subsampleauflösung 
wie die Phasenverschiebung aber sollte doch gehen. Die 
KorrelationsFUNKTION läßt sich schnell über drei FFTs rechnen:

Correlation = IFFT(FFT(x)*conjugiertkomplex(FFT(y))).

FFTs für den dsPic wirds schon fertig geben.

Ist doch nen schönes Projekt, ich hatte mal ne 3D Maus erstanden, die 
wurde durch 3 Mics im Raum geortet, hat solala gefunzt.

Cheers
Detlef

[Purzel H. (hacky)]

Eine fouriertrafo fuer eine einzelne Frequenz ist overkill. Da genuegt 
es mit dem sin und dem cos dieser frequenz zu multiplizieren.

[Tony P. (micropower)]

Ah Ha schrieb:
> Eine fouriertrafo fuer eine einzelne Frequenz ist overkill.

Richtig, deshalb möchte ich mein Glück mit dem Goertzel-Algorithmus 
probieren.

"imag" und "real" für eine Frequenz, sind mit dem Goertzel-Algorithmus 
schnell berechnet. Dann bekomme ich pro Kanal diese zwei Werte.

Jetzt, um die Richtung der Signale ausfindig zu machen, benötige ich 
theoretisch mindestens zwei Kanäle, aber eigentlich habe ich drei 
Kanäle.

Meine Frage ist; muss ich diese atan2 wirklich dreimal berechnen (einmal 
pro Kanal), oder kann ich die Imag- und Real-Werte vorher irgendwie 
"zusammenfassen" und mit eins oder zwei atan2-Berechnungen auskommen?