Goertzel Phasenbestimmung

Hallo,

ich versuche gerade verzweifelt den Görtzel Algorithmus auf einem 
Microcrontoller zu implemetieren. Habe aber noch Verständnisfragen.

Ich möchte nur die Phase von zwei Sinus Signalen auslesen. Beide Signale 
liegen bei etwa 25 kHz. Das Original Signal ist ein Sinus mit einer 
Amplitude von ca. 1,25 Volt. Das verschobene Signal ein Sinus der etwa 
um 40 bis 60 Grad verschoben ist und eine Amplitude von ca. 0,8 Volt 
hat. Und ich taste mit 120 kHz ab. Ich gebe immer 100 Abtastwerte auf 
den Algorithmus. Ich möchte die Differenz der Phasenlagen berechnen. 
Jedoch ändert sich diese ständig, obwohl die Signale immer die gleiche 
Phasenlage zueinander haben.

Was muss ich beachten bezüglich Abtastwerte oder besser gesagt 
Blockgröße?

Danke für eventuelle Tipps!

Gruß Marcus

Hallo,

@Detlef: danke für den Hinweis. Werde dies berücksichtigen. Hängt dies 
dann nicht auch an der benötigten Bandbreite. Was wäre wenn ich 120 Hz 
Bandbreite annehme?

@Alex: danke für den Hinweis. Sollte ich eventuell sogar das Signal mit 
diversen Fensterfunktionen (Hamming usw.) multiplizieren oder ist der 
Algorithmus auch mit einem normales Rechteckfenster zu gebrauchen?

Danke für die Hinweise.

Gruß Marcus

Alles klar danke. Werde es gleich mal mit Matlab simulieren.
Das mit der Fensterfunktion war klar. Meinte nur ob ich eventuell das 
Signal an den Rändern nicht so "hart" (Rechteckfenster) fenstern soll, 
sondern durch Hamming usw.

Gruß Marcus

Hallo Detlef,

ich hätte da noch eine Frage:

Ich erzeuge in Matlab einen Sinus mit der Frequenz 25 kHz und
nehme von diesem alle 4*10^-6 sec (entspricht einer Samplefrequenz von 
250 kHz) eine Stützstelle. Dann muss ich schauen, dass ich ein 
ganzzahliges Vielfaches von Perioden im "Fenster" habe. Dies wären doch 
dann n*T = n* 40*10^-6. Bandbreite = fs/N. Wie groß sollte die 
Bandbreite maximal sein? Und wenn ich an der Samplefrequenz nichts 
ändern möchte, muss ich ich die Blockgröße N dann auf z.B. auf 1250 
erhöhen? Dann hätte ich eine Bandbreite von 200 Hz. Und in N=1250 würden 
125*T hineinpassen.

Gruß Marcus

Nur so'n Kommentar, wo ich grade mitlese: Unter dem Stichwort "Sliding 
Goertzel" gibt's einige Erkenntnisse zu Phase und beliebigen Frequenzen 
(müssen nicht ins "Fenster" passen. Zeitliche Änderungen und laufende 
Ermittlung der Phase ist also ansich kein Problem.

Hallo zusammen,
hallo Detlef,

ich habe den Görtzel nun in Matlab realisieren können. Nun geht es 
darum, diesen Code auch im Controller zu implementieren.

Bis jetzt habe ich folgenden Code:

1

        double Fs = 229840;

2

  double k;

3

  double coeff;

4

  double cos_coeff;

5

  double realteil;

6

  double imagteil;

7

  double phi_rad;

8

  double phi_grad;

9

  double s0 = 0.0;

10

  double s1 = 0.0;

11

  double s2 = 0.0;

12

        int ind;

13

14

  k = round(f/Fs*N)+1; 

15

  coeff = 2*M_PI*(1/N)*k;

16

  cos_coeff = cos(coeff) * 2;

17

18

  x[N+1]=0;

19

  for (ind=1; ind==N+1; ind++) 

20

          {

21

      s0 = x[ind] + cos_coeff*s1 - s2;

22

          s2 = s1;

23

          s1 = s0;

24

      }

25

  realteil = s1 - s2*cos_coeff/2;

26

  imagteil = s2*sin(coeff);

27

  phi_rad= atan(imagteil/realteil);

28

29

  if (realteil < 0)

30

    {

31

    phi_rad = phi_rad - (M_PI/2);

32

    }

33

    else 

34

    {

35

    phi_rad = phi_rad + (M_PI/2);   

36

    }

37

38

  phi_grad = 180*phi_rad/M_PI;

39

40

41

  if (phi_grad > 180)

42

      {phi_grad = phi_grad - 360;}

43

        if (phi_grad < -180)

44

      {phi_grad = phi_grad +360;}

45

46

47

  return phi_grad;

Jedoch erhalte ich keine gescheiten Zahlen. Liegt wohl an den 
Nicht-Integer Berechnungen. Aber ich dachte, dass ich dies mit double 
Variablen richitg mache.
Kann mir da jemand auf die Sprünge helfen?
Danke.
Gruß Marcus

Hi Alex, SFT.., Detlef,

sorry, dass ich mich erst jetzt melde, hatte letzte Woche keine Zeit 
mehr.

Ja ich habe diesen Code in Matlab so laufen und er tut auch was er soll. 
Nur im Controller gab er mir die Phasen entweder gar nicht aus oder 
falsch. Ich habe mittlerweile den Code im Controller soweit am Laufen, 
dass er nun auch die Phase ausgibt. Nur habe ich bei der Berechnnung 
irgendwelche Schwankungen, so dass er immer zwischen 10 bis 20 Grad 
schwankt, aber stabil. Dieselbe Hardware nur mit Matlab über UART an Com 
Schnittstelle angesteuert und ausgelesen gibt mir exaktere Werte. D.h. 
es muss an der Berechnung im Controller liegen. Deswegen habe ich die 
float Zahlen in meiner Berechnung unter verdacht.
Hat da eventuell schon jemand den Görtzel auf dem Controller am Laufen 
mit Integer Zahlen (Code in C)?
Den Code habe ich gerade auf dem anderen Rechner und kann diesen 
deswegen jetzt nicht anhängen, aber ich werde dies morgen früh noch 
nachholen.


Gruß Marcus

Hallo Detlef,

danke für die tolle Hilfe.
Die Lösung meines Problems haben wir nun herausgefunden. Der 
µ-Controller
hatte Probleme bei der Abtastung. Zuerst war der Kerntakt zu langsam 
eingestellt und dann der Timer für die Abtastung durch den ADC.

Die Schwankungen sind soweit weg, da ich die Signale nun auch gescheit 
abtaste.

@alle anderen:
Wer diesen Thread zufällig liest und auch mit dem AdµC 7128/ 7129 
arbeitet und Hilfe braucht. Ich habe den C-Code nun fertig. Einfach 
anfragen.

Gruß Marcus

Hallo,

welchen unterschied würde es machen, wenn ich keine Blöcke durch den 
Filter schicke, sondern jedes Sample einzeln berechne.

Gruß

Also irgendwie ist das Ergebniss sehr unbefriedigend.
Ich möchte ja quasi den maximalen Amplitudenwert einer Sinusschwingung 
haben und die Phasenverschiebung dieser Schwingung gegenüber einem 
absoluten Wert.
Ob nach einem oder 8 Samples ist im Grunde egal, nur mehr oder weniger 
schnell sollte es sein..

Gruß QB

Ok,

dann ist klar, dass das Ergebniss unbefriedigend ist.
Abtasten tue ich mit 8,3Mhz. Der abzutastende Sinus hat 1Mhz.
Wie kann ich die Anzahl der benötigten Samples den bestimmen?
Ich benötige nur die 1Mhz alles andere kan wegfallen bzw. werde ich 
später versuchen dieses mit einem Vorgeschalteten bandpass noch 
rauszufiltern.

Gruß

Ok,

ich glaube so langsahm blicke ich dahinter.
mit der Anzazhl der Samples gebe ich im Brinzip meine Bandbreite vor. 
Die Frage ist, worauf bezieht sich diese?
Bei 8,3Mhz Abtastrate, und 100 Samples, wäre meine bandbreit somit 
800Hz?
Also 0 - 800Hz, oder 1Mhz +- 400Hz? Irgendwie leicht verwirrend.

Spielt es eine Rolle an welche Stelle des Signals ich Abtaste?
Ich meine, muss quasi synchron abgetastet werden, oder reicht es das die 
Anzahl der Samples pro periode gleich ist?

Ich meinte nur, wenn ich nicht genau mit einem Vielfachen von F0 
abtaste, verschieben sich meine Samples ja von periode zu periode.

Vielleicht habe ich auch kein Problem und ich verstehe auch was du mir 
sagst... also:

Sagen wir ich Taste immer mit einem Mehrfachen meiner Messfrequenz ab.
meine Abtastwerte sind somit bei jeder Periode z.B. bei:
1/4Pi, 1/2Pi, pi, 3/4 Pi.... es wiedreholt sich bei jeder periode.
Nun Taste ich um ein wenig mehr, als ein Mehrfaches der Abtastfrequnz 
ab:
(1/4+0.1)Pi, (1/2+0.2)Pi, pi+0.3, 3/4(0.4)Pi...
bei der nächsten Periode wäre dann mein erster abtastwert nicht mehr bei 
(1/4+0.1)Pi, sondern bei 1/4+0.5)Pi....
Für eine normale Abtastung spielt es keine Rolle, ich bin mir nur nicht 
sicher ob ich den Goerzel richtig verstanden habe.

Wie kann ich mir das berechnete Spektum anschauen? Wären das Betrag und 
Phase nach jedem Schleifendurchlauf?

Ich hab immernoch Probleme eine Plausible Amplidue und Phase zu 
bestimmen, die Werte springen stark. Hier mein Code:

#define F0    1000000    // Abgetastete Freq
#define FS    5000000    // Sampelfreq
#define T0    ((float)1/F0)
#define TS    ((float)1/FS)

void init_goerzel(_u32 n)
{
  _u32 k;
  float w;

  N = n*(T0/TS);
  B= FS/N;
  k = 0.5 + N*F0/FS;
  w = ((2*PI)/N)*k;
  cosine = cos(w) ;
  cos_coeff = 2*cosine;
  sine = sin(w);
}

void goerzel(_u16* x)
{
  float realteil;
  float imagteil;
  float phi_rad;
  float s0;
  float s1 = 0;    // bei jedem Block neu initialisieren
  float s2 = 0;    // bei jedem Block neu initialisieren
  int i =0;
  static int j=0;


  for (i=0;i<N;i++)
  {
    s0 = x[i] + cos_coeff*s1 - s2;
      s2 = s1;
      s1 = s0;
  }

    realteil = s1 - s2*cosine;
    imagteil = s2*sine;

  amplitude = sqrt(pow(realteil,2) + pow(imagteil,2));
    phi_rad= atan(imagteil/realteil);

    phi_grad = (phi_rad*180)/PI;
}

Derzeit ist die Blockgröße 15 also 3 Perioden.

Im ergebniss habe ich aber irgendwie nur schwachsinn...

Stimmt, habe ich schon gefunden den fehler, die Phase scheint nun zu 
stimmen, leider noch ein wenig verrauscht, ich muss noch ein wenig mit 
der Bandbreite Spielen. Leider habe ich nicht allzuviel Zeit um eine 
große Anzahl an Samples zu bestimmen am liebsten wären mit Amplitude und 
Phase schon nach einer Periode, ansonsten muss ich über eine 
Spitzenwertdetektion in Hardware nachdenken, was mir die möglichkeit 
rauben würde das Signal noch Digital zu Filtern... mal sehen.
Mit der Phase habe ich noch Probleme, aber ich glaube meine Abtastblöcke 
sind nicht Synchron. Ich werde heute versuchen dieses Problem zu lösen.

Ich hab mal ein Paar Messergebnisse angehängt.
Also die Schwankungen in der Amplitude mögen vorhanden sein, aber der 
Rest erscheint mir sehr verwirrend.
Ich nehme immer 32 Werte auf und beginne die Aufnahme dan von neuem. Ich 
verstehe nicht, warum die Phase so wandert.
Die sichbaren grenzen im Real- und Imaginärteil sind in der Berechnung 
nicht vorhanden. Mein Zahlenbereich hat für die Ausgabe nur nicht 
ausgereicht.

So, hab noch ein wenig was ausgebügelt, aber meine Phase ist ein 
Sägezahn. Einer ne Idee?

Könnte sein, dass meine Abtastfrequenz nicht genau ein vielfaches meines 
Signals ist oder?

Nun es scheint tatsächlich zu zu sein.
Kann ich den Filter auch als eine art Bandpass sehen? Oder brauche ich 
zusätzliche Filter um weitere Störungen auszugleichen?

Hm, also irgendwo muss da noch ein Wurm drin sein.
Auch wenn ich die Frequenz nachjustiere, komme ich nie genau auf den 
Wert, dass die Phase stehen bleibt. Mache ich eine "normale" FFT ist 
meine Phase bei 10Mhz ca. bei -3°, was auch der richtige wert sein 
Sollte. Mit dem Goerzel wandert sie aber.
Ich habe überlegt für diese stetige wanderung eine korrektufaktor 
einzubauen, gefallen tut mir dies aber nicht.

Kann es an einer "noch" fehlenden Fensterung liegen?

Mr. Qb schrieb:
> Kann es an einer "noch" fehlenden Fensterung liegen?

Nein, wenn du die Bedingungen bezüglich Blocklänge, Abtastung, 
Bandbreite ... einhältst ist es ja schon ein Rechteckfenster.

So richtig verstehe ich die Phasenberechnung vom Goerzel nicht.
Warum muss die Abtastfrequenz ein vielfaches der Samplefrequenz sein? 
Auch wenn ich asynchron abtaste, bleibt die Phase doch immer gleich. Es 
macht für mich keinen Sinn.

Hm ich verstehe es doch... lößt leider nicht mein Problem...

MrPPI schrieb:
> Hat jemand diesbezüglich eine Idee?

Beitrag "Re: Goertzel Phasenbestimmung"

Noch ein Nachtrag

Ohne deinen Code zu verstehen, wenn dein Abtastblock z.B. 50 Werte 
enthät, muß der 51(!) Wert Null enthalten und das Filtervon 0 bis 51 
laufen.

Gruß Joe

Dann sollte es laufen...