mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Goertzel-Algorithmus und Signed Fractional Format


Autor: Tony P. (micropower)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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!

Autor: Tony P. (micropower)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Detlef _a (detlef_a)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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 ;

}

Autor: Tony P. (micropower)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: Detlef _a (detlef_a)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
>>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

Autor: Zwölf Mal Acht (hacky)
Datum:

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

Autor: Tony P. (micropower)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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?

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.