Görtzel Algorithmus, falsche Werte bei DC Offset?

[Friedrich (Gast)]

Hallo an Fourier Experten!

ich benutze folgenden Code zum Testen des Görtzel Algorithmus:
https://github.com/OmaymaS/DTMF-Detection-Goertzel-Algorithm-

Das funktioniert soweit recht gut, allerdings wird mit steigerung des DC 
offsets das Ergebnis immer schlechter, legt man im Extremfall Aref an 
den ADC Eingang, das heißt alle Samplewerte sind das maximum (1023), ist 
das Ergebnis ganz furchtbar.

Das ist das Ergebnis der einzelnen Spektralwerte, der zu errechneten 
Frequenzen, bei ADCin = Aref
697Hz - 78
770Hz - 2828
852Hz - 3343
941Hz - 2455
1209Hz - 85
1336Hz - 1068
1477Hz - 597
1633Hz - 666

Eigentlich sollte der DC Offset keinen Einfluss auf die einzelnen 
Spektrallinien, die der Görtzel Algorithmus errechnet haben.

Hier ein Beispielergebnis bei 200mV Offset, 400mVpp 852Hz Sinus (Aref = 
1.1V)

697Hz - 38
770Hz - 57
852Hz - 817
941Hz - 25
1209Hz - 9
1336Hz - 3
1477Hz - 1
1633Hz - 2


Habe ich etwas übersehen?

Danke

[Georg G. (df2au)]

Friedrich schrieb:
> Habe ich etwas übersehen?

Welchen Wert liefert dein ADC, wenn der Eingang auf ARef oder darüber 
liegt? Gibt es da noch eine Abhängigkeit von der überlagerten DTMF 
Frequenz? Was sagt dein ADC Eingang dazu, dass du ihn in den positiven 
Halbwellen überfährst? Wie soll da noch ein vernünftiges Ergebnis 
kommen?

Der Wert der Eingangsspannung darf zu keiner Zeit kleiner als 0 oder 
größer als ARef sein.

[Friedrich (Gast)]

Hallo Georg,

Georg G. schrieb:
> Friedrich schrieb:
>> Habe ich etwas übersehen?
>
> Welchen Wert liefert dein ADC, wenn der Eingang auf ARef oder darüber
> liegt? Gibt es da noch eine Abhängigkeit von der überlagerten DTMF
> Frequenz?

Der liefert 1023 bei ADCin = Aref, bzw. aufgrund des Shifts um 2 
Stellen, 255 am Eingang des Görtzel Algorithmus. Nein da gibt es keine 
Abhängigkeiten der DTMF, alle ADC Samples entsprechen dann 255.

> Was sagt dein ADC Eingang dazu, dass du ihn in den positiven
> Halbwellen überfährst? Wie soll da noch ein vernünftiges Ergebnis
> kommen?

Wo überfahre ich ihn mit positiven Halbwellen? Ich habe zum Test den 
ADCin auf Aref gelegt, und ihn nicht mit zu hoher Spannung oder dgl. 
betrieben.

> Der Wert der Eingangsspannung darf zu keiner Zeit kleiner als 0 oder
> größer als ARef sein.

Richtig, deshalb auch ADCin = Aref beim Test

[Audiomann (Gast)]

Ich kenne den Görtzel nicht genau, weiß aber aus Erfahrung dass FFT auch 
auf offset reagiert. Der Offset ist ein Rechteck mit hohen Frequenzen, 
die von der FFT mitgenommen werden, da sie wie auch die künstlichen 
Frequenzen durch das leakage beim Fenstern nicht völlig weggenommen 
werden.

Ich nehme an, das ist beim Görtzel genau so?

Frage: Kann man nicht einfach die Werte, die man in den Görtzel wirft, 
zuvor durch Mittelung und Abziehen vom Offset befreien?

[Wolfgang (Gast)]

Audiomann schrieb:
> Der Offset ist ein Rechteck mit hohen Frequenzen,
> die von der FFT mitgenommen werden, da sie wie auch die künstlichen
> Frequenzen durch das leakage beim Fenstern nicht völlig weggenommen
> werden.

Der Offset ist eine DC-Komponente, die nach der FFT bei 0Hz auftauchtt.
Was meinst du mit "hohen Frequenzen"?

[Dergute W. (derguteweka)]

Moin,

Wolfgang schrieb:
> Der Offset ist eine DC-Komponente, die nach der FFT bei 0Hz auftauchtt.
> Was meinst du mit "hohen Frequenzen"?

Ja, in der Theorie, aufm Stueckchen Karopapier. Aber in echten CPUs 
koennen halt mal Ueberlauefe auftauchen. Und genau nach sowas siehts mir 
aus. Der ganze Code ist nicht grad' vertrauenserweckend. Allein schon 
solche Konstrukte:

int N = 96;              // block size
volatile int samples[96]; // buffer to store N samples

staerken mein Misstrauen gegen den Rest ganz erheblich.
Wenn man da jetzt mal genau anfaengt, bei jeder Multiplikation/Addition 
zu stoebern, was da schief gehen kann, dann wird man sicher fuendig 
werden.

Gruss
WK

[Friedrich (Gast)]

Audiomann schrieb:
> Frage: Kann man nicht einfach die Werte, die man in den Görtzel wirft,
> zuvor durch Mittelung und Abziehen vom Offset befreien?

Hallo, ja genau das habe ich mittlerweile gemacht, und damit sind die 
Ergebnisse auch sehr gut.

Trotzdem würde mich interessieren, ob es am Algorithmus an sich liegt, 
oder ob da ein Fehler vorhanden ist.

Dergute W. schrieb:
> Moin,
>
> Wolfgang schrieb:
>> Der Offset ist eine DC-Komponente, die nach der FFT bei 0Hz auftauchtt.
>> Was meinst du mit "hohen Frequenzen"?
>
> Ja, in der Theorie, aufm Stueckchen Karopapier. Aber in echten CPUs
> koennen halt mal Ueberlauefe auftauchen. Und genau nach sowas siehts mir
> aus. Der ganze Code ist nicht grad' vertrauenserweckend. Allein schon
> solche Konstrukte:int N = 96;              // block size
> volatile int samples[96]; // buffer to store N samplesstaerken mein
> Misstrauen gegen den Rest ganz erheblich.

Ja das ist leider etwas unglücklich geschrieben, der Rest macht mir (als 
Hobbyprogrammierer) aber einen recht soliden Eindruck.

> Wenn man da jetzt mal genau anfaengt, bei jeder Multiplikation/Addition
> zu stoebern, was da schief gehen kann, dann wird man sicher fuendig
> werden.
>
> Gruss
> WK

[c-hater (Gast)]

Dergute W. schrieb:

> Wolfgang schrieb:
>> Der Offset ist eine DC-Komponente, die nach der FFT bei 0Hz auftauchtt.
>> Was meinst du mit "hohen Frequenzen"?
>
> Ja, in der Theorie, aufm Stueckchen Karopapier. Aber in echten CPUs
> koennen halt mal Ueberlauefe auftauchen.

Das ist der Punkt: man muss die Algorithmen so auslegen, dass auch im 
worst case keine Überläufe auftreten, sonst kommt Müll raus. Der worst 
case für einen Goertzel ist übrigens nicht DC, sondern seine 
Resonanzfrequenz...

Sprich: Wenn's schon bei DC knallt, ist der Goertzel definitiv falsch 
implementiert! Ein korrekt implementierter Goertzel liefert bei DC 
natürlich 0 (oder sehr nahe dran) und nix anderes.

Wobei der durch DC verursachte Sprung schon zu berücksichtigen ist. Der 
liefert immer auch Frequenzbestandteile, die die Resonanzfrequenz des 
Görtzel treffen. Die darin enthaltene Energie findet sich natürlich am 
Ende im Goertzel-Ergebnis wieder. Man kann das durch geeignete 
Initialisierung der Goertzel-History oder durch eine vorgeschaltete 
Fensterfunktion aber vollständig oder zumindest weitgehend unterdrücken.

Wie schon gesagt: DC ist unproblematisch, selbst von einem nicht 
kompensierten Sprung durch DC kommt logischerweise umso weniger am 
Ausgang an, je schmalbandiger der Goertzel ist.

Überläufe sind bei der Resonanzfrequenz am wahrscheinlichsten! Die 
Implementierung muss so ausgelegt sein, dass auch bei Resonanz noch kein 
Überlauf eintritt.

[Dergute W. (derguteweka)]

Moin,

Friedrich schrieb:
> der Rest macht mir (als
> Hobbyprogrammierer) aber einen recht soliden Eindruck.

Q = (sample[i]) + ((coeff * Q_prev) >> 14) - (Q_prev2); // >>14 used as the coeff was used in Q15 format

Also z.B. in obiger Zeile hab' ich ein leichtes Ziehen im Bauch:
coeff * Q_prev - das sind beides (16bit?) Integer, also kommt als 
Produkt auch ein Integer raus. Und das wird dann 14 bit nach rechts 
geschoben...
Neeneeneee...

Gruss
WK

[c-hater (Gast)]

c-hater schrieb:

> Überläufe sind bei der Resonanzfrequenz am wahrscheinlichsten! Die
> Implementierung muss so ausgelegt sein, dass auch bei Resonanz noch kein
> Überlauf eintritt.

Da fehlte noch der Halbsatz:

...dann funktionieren sie auch bei allen Frequenzen 0..<1/2fs 
zuverlässig, also insbesondere einschließlich DC.

[Friedrich (Gast)]

Dergute W. schrieb:
> Moin,
>
> Friedrich schrieb:
>> der Rest macht mir (als
>> Hobbyprogrammierer) aber einen recht soliden Eindruck.
> Q = (sample[i]) + ((coeff * Q_prev) >> 14) - (Q_prev2); // >>14 used as
> the coeff was used in Q15 format
> Also z.B. in obiger Zeile hab' ich ein leichtes Ziehen im Bauch:
> coeff * Q_prev - das sind beides (16bit?) Integer, also kommt als
> Produkt auch ein Integer raus. Und das wird dann 14 bit nach rechts
> geschoben...
> Neeneeneee...

Naja, soweit ich verstanden habe, kommen die 14 Bit kommen daher, dass 
dort keine Fließkommaarithemtik verwendet wird, er bläst ja auch die 
Koeffizienten um (1<<14) auf, deshalb muss er das auch wieder raus 
bringen.

c-hater schrieb:
> Dergute W. schrieb:
>
> Das ist der Punkt: man muss die Algorithmen so auslegen, dass auch im
> worst case keine Überläufe auftreten, sonst kommt Müll raus. Der worst
> case für einen Goertzel ist übrigens nicht DC, sondern seine
> Resonanzfrequenz...
>
> Sprich: Wenn's schon bei DC knallt, ist der Goertzel definitiv falsch
> implementiert! Ein korrekt implementierter Goertzel liefert bei DC
> natürlich 0 (oder sehr nahe dran) und nix anderes.

Hallo c-hater,
siehst du in der Implementierung einen Fehler?
Für mich sieht es schlüssig aus, ohne Fehler in der Implementierung des 
Görtzel Algroithmus.

goertzel (int sample[], long int coeff, int N)
//---------------------------------------------------------------//     
{
//initialize variables to be used in the function
  int Q, Q_prev, Q_prev2, i;
  long prod1, prod2, prod3, power;

  Q_prev = 0;      //set delay element1 Q_prev as zero
  Q_prev2 = 0;      //set delay element2 Q_prev2 as zero
  power = 0;      //set power as zero

  for (i = 0; i < N; i++)  // loop N times and calculate Q, Q_prev, Q_prev2 at each iteration
    {
      Q = (sample[i]) + ((coeff * Q_prev) >> 14) - (Q_prev2);  // >>14 used as the coeff was used in Q15 format
      Q_prev2 = Q_prev;    // shuffle delay elements
      Q_prev = Q;
    }

  //calculate the three products used to calculate power
  prod1 = ((long) Q_prev * Q_prev);
  prod2 = ((long) Q_prev2 * Q_prev2);
  prod3 = ((long) Q_prev * coeff) >> 14;
  prod3 = (prod3 * Q_prev2);

  power = ((prod1 + prod2 - prod3)) >> 8;  //calculate power using the three products and scale the result down

  return power;
}

[Dergute W. (derguteweka)]

Moin,

Friedrich schrieb:
> Naja, soweit ich verstanden habe, kommen die 14 Bit kommen daher, dass
> dort keine Fließkommaarithemtik verwendet wird, er bläst ja auch die
> Koeffizienten um (1<<14) auf, deshalb muss er das auch wieder raus
> bringen.

Ja, klar - aufm Papier stimmt ja die Rechnung. Aber bei der Schreibweise 
in C kommt bei einer 16x16 bit Multiplikation auch nur ein 16 bit 
Produkt raus. Da kanns dann den Ueberlauf geben.

Gruss
WK

[Rainer V. (a_zip)]

Jedenfalls ist klar, dass der Gleichspannungsoffset bei der FFT "in" der 
0-Linie erscheint, wenn er den Eingang nicht in einen unzulässigen 
Bereich verschiebt. Vielleicht kannst du dein Signal ja mal über einen 
Kondensator einkoppeln und die FFT vergleichen...kann man übrigens für 
Untersuchungen auch sehr schön mit einem Simulationsprogramm anschauen!
Gruß Rainer

[c-hater (Gast)]

Friedrich schrieb:

> siehst du in der Implementierung einen Fehler?

Einen? Das wäre ja schon richtig gut für einen C-ler..

Es fängt doch schon in der Funktionsdeklaration an. Da wird N übergeben. 
Das aber ist ein sehr wesentliches Kriterium dafür, was an 
Verarbeitungsbitbreite letztlich nötig sein wird...

Aber geht der Code in irgendeiner Form darauf ein? Nein, wie von 
stupiden C-lern kaum anders zu erwarten, tut er das natürlich nicht.

Oder anders ausgedrückt: Diese Funktion ist bereits vom Entwurf her 
kompletter Schwachsinn.

Und da wundern sich die Leute, warum ich C so hassenswert finde...

Das verklebt den Leuten nur den Blick für's Wesentliche, durch 
Vortäuschung einer objektiv nicht vorhandenen Typsicherheit. Was nützt 
denn ein Typ, dessen Grenzen ich jederzeit ohne jede Warnung verletzen 
kann?

[Friedrich (Gast)]

Dergute W. schrieb:
> Moin,
>
> Friedrich schrieb:
>> Naja, soweit ich verstanden habe, kommen die 14 Bit kommen daher, dass
>> dort keine Fließkommaarithemtik verwendet wird, er bläst ja auch die
>> Koeffizienten um (1<<14) auf, deshalb muss er das auch wieder raus
>> bringen.
>
> Ja, klar - aufm Papier stimmt ja die Rechnung. Aber bei der Schreibweise
> in C kommt bei einer 16x16 bit Multiplikation auch nur ein 16 bit
> Produkt raus. Da kanns dann den Ueberlauf geben.
>
> Gruss
> WK

Hallo,
habe jetzt alles auf Float umgebaut, und komme leider zum gleichen 
ernüchternden Ergebnis. Ebenso mit dem Code aus diesem Beitrag: 
https://stackoverflow.com/questions/11579367/implementation-of-goertzel-algorithm-in-c
Auch hier ist der Wert bei DC Offset falsch. Rechne ich den DC Offset 
raus, funktioniert der Algorithmus.

Rainer V. schrieb:
> Jedenfalls ist klar, dass der Gleichspannungsoffset bei der FFT
> "in" der
> 0-Linie erscheint, wenn er den Eingang nicht in einen unzulässigen
> Bereich verschiebt. Vielleicht kannst du dein Signal ja mal über einen
> Kondensator einkoppeln und die FFT vergleichen...kann man übrigens für
> Untersuchungen auch sehr schön mit einem Simulationsprogramm anschauen!
> Gruß Rainer

Du meinst, dass ich von meinem Signal einmal eine vollständige FFT 
machen soll, damit ich auch den DC Beitrag explizit errechne?

c-hater schrieb:
> Friedrich schrieb:
>
>> siehst du in der Implementierung einen Fehler?
>
> Einen? Das wäre ja schon richtig gut für einen C-ler..
>
> Es fängt doch schon in der Funktionsdeklaration an. Da wird N übergeben.
> Das aber ist ein sehr wesentliches Kriterium dafür, was an
> Verarbeitungsbitbreite letztlich nötig sein wird...
>
> Aber geht der Code in irgendeiner Form darauf ein? Nein, wie von
> stupiden C-lern kaum anders zu erwarten, tut er das natürlich nicht.
>
> Oder anders ausgedrückt: Diese Funktion ist bereits vom Entwurf her
> kompletter Schwachsinn.
>
> Und da wundern sich die Leute, warum ich C so hassenswert finde...
>
> Das verklebt den Leuten nur den Blick für's Wesentliche, durch
> Vortäuschung einer objektiv nicht vorhandenen Typsicherheit. Was nützt
> denn ein Typ, dessen Grenzen ich jederzeit ohne jede Warnung verletzen
> kann?

Sorry c-hater, aber irgendwie erkenne ich darin keinen Ansatz der 
konstrukiven Lösungsfindung?
Du hast natürlich Recht, aber in dem spezifischen Fall von N = 96 ist 
alles gegeben um keinen Überlauf zu riskieren, also ist das kein Fehler 
betreffend des DC-Offset Problems.

[Wolfgang (Gast)]

Dergute W. schrieb:
> Ja, klar - aufm Papier stimmt ja die Rechnung. Aber bei der Schreibweise
> in C kommt bei einer 16x16 bit Multiplikation auch nur ein 16 bit
> Produkt raus.

Was passiert wohl, wenn der Compiler das "long int" vor coeff sieht?
Es wäre "ungewöhnlich", wenn er dann zu einer 16x16 Multiplikation 
greift.
Oder was meinst du genau?

[Dergute W. (derguteweka)]

Moin,

Ja, haste Recht. Weiter oben ist coeff nur "normaler" int, dass das 
nochmal als long int uebergeben wird, hab' ich uebersehen. Aber nach 
irgendsowas von dem Kaliber riecht's mir schon. Bin aber nicht motiviert 
genug, mich da durchzufieseln.

Gruss
WK

[El Ef (Gast)]

Ich denke dass liegt eher an den Testbedingungen die mir hier nicht ganz 
klar sind. Hast du mal probiert dir dein Signal am Rechner zu erzeugen 
mit und ohne offset und den Algorithmus damit bespaßen? Hast du dir 
deine Adc Samples mal aufgezeichnet und angeschaut ob das dass so 
aussieht wie erwartet?

[Rainer V. (a_zip)]

Friedrich schrieb:
> Du meinst, dass ich von meinem Signal einmal eine vollständige FFT
> machen soll, damit ich auch den DC Beitrag explizit errechne?

Ja. Und ich meine, du solltest das mit der C-Ankopplung mal 
probieren...wenn du schon nicht simulieren willst...(mir ist allerdings 
nicht klar, wie man den Görtzel simulieren könnte...)
Also genug zum Probieren, viel Spass und Gruß, Rainer

[Friedrich (Gast)]

Leider sind meine Mathcad Kentnisse recht eingestaubt, aber soweit 
sollte es reichen.

(Getestet habe ich die Funktion aus dem 2. Link von stackoverflow)

Auch mittels Mathcad erzeugt er einen Output bei, hier gegebenen 770Hz.

Alle Werte die Eingangsseitig verarbeitet werden sind "255".

Zum Vergleich eine DFT mit 16 Samples - hier kommt wie zu erwarten nur 
DC Anteil raus.

[Friedrich (Gast)]

Ich denke den Fehler gefunden zu haben:

k does not need to be an integer.  But if it isn't, there will
be aliasing from other frequencies whose periods are exact
submultiples on N, including DC (which isn't the case when
k is an integer).  If this difference could cause any problems,
windowing before the Goertzel sum might help minimize some
of the results of this difference.

Runde ich k auf eine Ganzzahl, dann sind die Ergebnisse im Bereich 
10^-12, also quasi 0.
Teste es Morgen noch in der Realität, nicht nur in der Simulation.

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
    int     k,i;
    float   floatnumSamples;
    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

    float   scalingFactor = numSamples / 2.0;

    floatnumSamples = (float) numSamples;
    k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
    omega = (2.0 * M_PI * k) / floatnumSamples;
    sine = sin(omega);
    cosine = cos(omega);
    coeff = 2.0 * cosine;
    q0=0;
    q1=0;
    q2=0;

    for(i=0; i<numSamples; i++)
    {
        q0 = coeff * q1 - q2 + data[i];
        q2 = q1;
        q1 = q0;
    }

    // calculate the real and imaginary results
    // scaling appropriately
    real = (q1 - q2 * cosine) / scalingFactor;
    imag = (q2 * sine) / scalingFactor;

    magnitude = sqrtf(real*real + imag*imag);
    return magnitude;
}

[Dergute W. (derguteweka)]

Moin,

Guter Punkt. Da haett' ich nicht als erstes drangedacht.

Gruss
WK

[Friedrich (Gast)]

Moin,

also das dürfte tatsächlich das Problem des DC Offsets lösen!

Hat aber die unangenehme Nebenwirkung, dass ich bei meinen derzeit 
gegebenen Daten, und der Annahme, dass k eine Ganzzahl sein muss:

N = 96
Samplingrate = 9615

 k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));

nur eine "Auflösung" von 100Hz habe (50Hz - 150Hz, ...)
Soll heißen dass z.B. zwischen den Frequenzen 852Hz und 941Hz keine 
Unterscheidung mehr stattfinden kann, da bei beiden Frequenzen k = 9 
ist.

Das hieße mein Verhältnis N / Samplingrate müsste erheblich größer 
werden, um ein Vernünftiges Ergebnis zu erzielen.

Die Besten Ergebnisse habe ich bislang erhalten, wenn ich den DC Offset 
berechnet und entfernt habe, bevor die Messwerte in den Görtzel 
Algorithmus gelangen (ohne k auf Ganzzahl zu runden).

Daher wäre die Frage ob es nicht so sinnvoller ist, oder aber zu 
Fenstern?

[Rainer V. (a_zip)]

Friedrich schrieb:
> Daher wäre die Frage ob es nicht so sinnvoller ist, oder aber zu
> Fenstern?

Und bei den "Fenstereien" sind wiederum jede Menge Einzelheiten zu 
berücksichtigen! Das ist nicht einfach. Und an 20 Messwerte mal locker 
530 Nullen anhängen, könnte ein "schönes" Ergebnis zeigen...aber dann 
mußt du definitiv wissen, was du da gemacht hast. Ich bleibe bei meinem 
Vorschlag, erst einmal eine einfache FFT mit und ohne Koppelkondensator 
und hast du denn überhaupt so was wie ein Referenzsignal, bei dem du 
weißt, was deine Software rausschmeissen muß? Jetzt meinetwegen an 
Rundungsfehlern rumzuschrauben bringt doch gar nichts!
Gruß Rainer

[Dergute W. (derguteweka)]

Moin,

Friedrich schrieb:
> Die Besten Ergebnisse habe ich bislang erhalten, wenn ich den DC Offset
> berechnet und entfernt habe, bevor die Messwerte in den Görtzel
> Algorithmus gelangen (ohne k auf Ganzzahl zu runden).

Naja, dann mach das doch so...
Vielleicht weniger Rechenaufwand: Wenn du statt der Samples selber immer 
nur die Differenz zum vorherigen Sample "goertzelst"; das wirkt dann wie 
ein Differenzierer, also eine Art Hochpass. Kann aber sein, dass du dann 
zuviel Rauschen im Signal hast.

Gruss
WK