Görtzel Algorithmus, falsche Werte bei DC Offset?

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

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?

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"?

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

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.

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.

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.

1

goertzel (int sample[], long int coeff, int N)

2

//---------------------------------------------------------------//     

3

{

4

//initialize variables to be used in the function

5

  int Q, Q_prev, Q_prev2, i;

6

  long prod1, prod2, prod3, power;

7

8

  Q_prev = 0;      //set delay element1 Q_prev as zero

9

  Q_prev2 = 0;      //set delay element2 Q_prev2 as zero

10

  power = 0;      //set power as zero

11

12

  for (i = 0; i < N; i++)  // loop N times and calculate Q, Q_prev, Q_prev2 at each iteration

13

    {

14

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

15

      Q_prev2 = Q_prev;    // shuffle delay elements

16

      Q_prev = Q;

17

    }

18

19

  //calculate the three products used to calculate power

20

  prod1 = ((long) Q_prev * Q_prev);

21

  prod2 = ((long) Q_prev2 * Q_prev2);

22

  prod3 = ((long) Q_prev * coeff) >> 14;

23

  prod3 = (prod3 * Q_prev2);

24

25

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

26

27

  return power;

28

}

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?

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.

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?

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?

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.

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.

1

float goertzel_mag(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)

2

{

3

    int     k,i;

4

    float   floatnumSamples;

5

    float   omega,sine,cosine,coeff,q0,q1,q2,magnitude,real,imag;

6

7

    float   scalingFactor = numSamples / 2.0;

8

9

    floatnumSamples = (float) numSamples;

10

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

11

    omega = (2.0 * M_PI * k) / floatnumSamples;

12

    sine = sin(omega);

13

    cosine = cos(omega);

14

    coeff = 2.0 * cosine;

15

    q0=0;

16

    q1=0;

17

    q2=0;

18

19

    for(i=0; i<numSamples; i++)

20

    {

21

        q0 = coeff * q1 - q2 + data[i];

22

        q2 = q1;

23

        q1 = q0;

24

    }

25

26

    // calculate the real and imaginary results

27

    // scaling appropriately

28

    real = (q1 - q2 * cosine) / scalingFactor;

29

    imag = (q2 * sine) / scalingFactor;

30

31

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

32

    return magnitude;

33

}

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

1

 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?