Hallo, hat jemand eigentlich schon mal die Phase eines Sinus über Goertzel richtig berechnet? Die Amplitude bekomme ich genau hin doch die Phase ist immer Mist. Ich habe mich genau an diese Vorschrift gehalten. http://www.eetimes.com/design/embedded/4024443/The-Goertzel-Algorithm Gruß D.
Das wird Dir ohne weitere Infos niemand beantworten können. Warum meinst Du ist die Phase nicht OK? Wie stellt sich das Problem genau dar? Berechnest Du kontinuierlich, oder blockweise? Wie sieht Dein Code aus?
Ich berechne Blockweise. Nach der Berechnung stimmt der Betrag der Amplitude exakt mit dem Sinus überein. Die Phase hat jedoch einen konstanten Fehler in Abhängigkeit der Blocklänge und der Abtastzeit.
>Nach der Berechnung stimmt der Betrag der >Amplitude exakt mit dem Sinus überein. Die Phase hat jedoch einen >konstanten Fehler in Abhängigkeit der Blocklänge und der Abtastzeit. Was willst du mit der Phase? Soweit ich Goertzel und FFT verstanden habe kann man damit feststellen ob eine Frequenz in den Samples vorhanden ist. Man kann auch feststellen mit welchem Energieinhalt. Man kann nicht feststellen zu welchem Zeitpunkt sie in den Samples vorhanden ist. Sie kann in der ersten Hälfte der Samples mit sehr hohem Pegel vorhanden sein und in der zweiten gar nicht. Wenn sie über alle Samples mit halber Amplitude vorhanden ist bekommst du das gleiche Ergebnis. Eine Phase kannst du damit nicht feststellen. In Bezug auch was auch?
Natürlich kann man mit FFT, Goertzel oder SDFT Algorithmen die Phase bestimmen. @Detlef: Ist der Fehler zufällig ein Vielfaches von "pi"? <der Holger mit Grossbuchstaben>
Holger schrieb: > @Detlef: Ist der Fehler zufällig ein Vielfaches von "pi"? Ja, danke! Die atan-Berechnung benötigt noch die Korrektur um pi. Das war aber nicht der Fehler. Der korrekte Amplitudenwert ist gegenüber den korrekten Phasenwert genau um ZWEI Abtastschritte verschoben.
Wieviel Werte shcibest Du denn rein und welche Frequenz interessiert dich?
Hi ich verwende den G. zum Berechnen von Phase und Amplitude und hab keine Probleme. Es geistern aber im Netz aber einige falsche Implementierungen herum. Häufig hakt es an folgenden Umständen: - die Speicher-Zellen müssen vor dem Starten mit 0 initialisiert sein. - im letzten Durchlauf muss dem Eingangssignal eine 0 angehängt werden. Um den G. einer Sequenz mit N Samples zu berechnen benötigt man also N+1 Durchläufe des Filters. Grüße, Alex
@Alex Danke! Das war es. Die Verschiebung von N-2 deutete irgendwie schon darauf hin. Nun läuft es perfekt. Detlef
Dank der kompetenten Hilfe hier im Forum nun ein Goertzelalgorithmus der fehlerfrei arbeitet. Für alle Wiederholungstäter - siehe Anhang ;-) Detlef
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 Marcus, deine Filterbandbreite ist zu groß. Du musst entweder die Abtastfrequenz kleiner machen oder die Blockgröße erhöhen. Anbei mal zwei Varianten. 1. 120 kHz Abtastfrequenz und Blockgröße 600 2. 20 kHz Abtastfrequenz und Blockgröße 100 Der Vorteil von Görtzel liegt ja gerade in der Unterabtastung.
Der Goertzel (wie auch der FFT) Algorithmus liefert nur korrekte Werte für ganzzahlige Frequenzen (1/2/3/4... Schwingungen pro Block (bzw. Fenster)). Falls die Frequenz nicht ganzzahlig im Fenster ist -> leakage effect ( http://de.wikipedia.org/wiki/Leck-Effekt ). Gruß, Alex
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
Marcus Moser schrieb: > Was wäre wenn ich 120 Hz > Bandbreite annehme? Das kannst du natürlich machen. Bei einer Abtastfrequenz von 20 kHz kommst du bei B=120 aber nicht auf ein ganzzahlige Blockgröße. Bei B=125 erhältst du jedoch für die Blockgröße N=160. Außerdem muss, wie Niedes schon sagte, immer eine ganzzahlige Periode der Schwingungen in den Block passen. Das ist bei N=160 auch erfüllt. Das „fenstern“ kannst du dir damit sparen, es ist ja somit ein Rechteckfenster.
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
Marcus Moser schrieb: > Meinte nur ob ich eventuell das > Signal an den Rändern nicht so "hart" (Rechteckfenster) fenstern soll Wenn du dafür sorgst, dass immer eine volle Periode oder ein Vielfaches davon in den Block passt (Phasenlage egal), dann musst du nicht fenstern. Wenn du nun virtuell die Blöcke ins positive und ins negative aneinander hängen würdest, entsteht an den Blockgrenzen keine Unstetigkeit. Damit ist existiert das Signal von Minus Unendlich bis Plus Unendlich. Somit tritt kein Leakage auf.
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
@Marcus Deine Überlegung ist richtig. Du musst nur die Bedingung n*T = N*Ts einfach einhalten. Ist wie in deinem Fall n=125 dann ist N=1250. Es geht natürlich auch mit n=1 und damit mit N=10. Allerdings wird die Bandbreite katastrophal. Praktisch würde ich zunächst meine maximale Bandbreite ermitteln und dann das zugehörige n bzw. N ermitteln.
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
Hallo! Also im Görtzel-Teil sehe ich auf die Schnelle keinen Fehler. Was meinst du mit "keine gescheiten Zahlen"? Hast du die Berechnung in Matlab genau so gemacht? Alex
Hi,
arbeitest Du mit MathCad? Das Ablaufdiagramm in den ersten Posts sieht
zumindest danach aus.
Wenn es MathCad ist, könnte es vielleicht am
atan (imag / real)
liegen, wenn die Phase "komisch" ist, da sich dieses "atan(...)" immer
auf den 1. Quadranten des komplexen Ebene bezieht. Wenn die komplexe
Zahl dort nicht liegt, produziert der "atan(...)" halt Käse.
Probier die Phase mal anders zu berechnen. Ich glaube, in MathCad gibt's
noch den Befehl "angle(...)" (oder so). Dabei wird automaisch der
richtige Quadrant bestimmt.
Dann klappt's auch mit der Phase ;-)
Marcus Moser schrieb: > Jedoch erhalte ich keine gescheiten Zahlen. Der Code sieht korrekt aus. Was ist dein Fehler? Gebe mal ein Bsp. an.
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
Ich habe gerade mal mit unterschiedlichen Rundungsvarianten gerechnet. 3 Nachkommastellen 0.1 % Phasenfehler 2 Nachkommastellen 1 % Phasenfehler 1 Nachkommastelle 4.5 % Phasenfehler 0 Nachkommastellen 18 % Phasenfehler
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
Mr. Qb schrieb: > Also irgendwie ist das Ergebniss sehr unbefriedigend. Dann hast du noch einen Fehler in deiner Rechnung. Der Algorithmus rechnet exakt. Mr. Qb schrieb: > Ob nach einem oder 8 Samples ist im Grunde egal, nur mehr oder weniger > schnell sollte es sein.. Das ist ein Widerspruch. Der Algorithmus benötigt je nach Abtastfrequenz und Filterbandbreite eine exakt darauf abgestimmte Anzahl von Sampels, irgendeine Anzahl zu nehmen produziert mit Sicherheit keine richtige Lösung.
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.
Du musst zunächst dafür sorgen, dass in einem Datenblock immer ein n-faches einer Periode ist. Hältst du diese Bedingung nicht ein, so springt deine Phase von Block zu Block. Weiterhin benötigst du eine entsprechende Bandbreite. Ein Goertzelfilter ist ja nur eine vereinfachte FFT mit genau einer Frequenz. Die Bandbreite bestimmt sozusagen die Filterbandbreite. Um es in Formeln zu pressen: f0 Frequenz welche detektiert werden soll, der Kehrwert T0 fs Sampelfrequenz, der Kehrwert TS n Anzahl der zu erfassenden Perioden Blocklänge: N = n * (T0/TS) Bandbreite: B = fs/N
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?
1. Die Stelle spielt eine Rolle, z.B. immer im Nulldurchgang wird natürlich nichts. 2. Grundvoraussetzung sind natürlich äquidistante Abtastschritte.
Ich meinte nur, wenn ich nicht genau mit einem Vielfachen von F0 abtaste, verschieben sich meine Samples ja von periode zu periode.
Vielleicht verstehe ich ja dein Problem nicht. Ein Block ist immer so aufgebaut, dass die Blöcke aneinander gehangen werden können und sich damit ein stetiges Signal ergeben würde, d.h. auf den letzten Sampelwert folgt wieder der erste Wert. Ansonsten bekommst du den Leakage effect. Ist diese Bedingung eingehalten, berechnet Goertzel exakt Amplitude und Phasenlage dieses Blockes. Sind im nächsten Block die Sampels verschoben, bekommst du selbstverständlich eine andere Phasenlage für diesen Block.
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.
Mr. Qb schrieb: > Für eine normale Abtastung spielt es keine Rolle, ich bin mir nur nicht > sicher ob ich den Goerzel richtig verstanden habe. Tatsächlich spielt es für die FFT auch eine Rolle. Mache dazu das einfache Gedankenmodell. Die abgetastete Funktion ist ein reiner Sinus also eine ungerade Funktion. Dann besitzt das Spektrum nur einen Imaginärteil. Ist die abgetastete Funktion ein reiner Cosinus, also eine gerade Funktion, dann besitzt das Spektrum nur einen Realteil. Verschiebt sich also dein Abtastwert leicht, ändert sich das Spektrum kontinuierlich zwischen Imaginärteil und Realteil. Wenn du aus diesen beiden Werten nur den Betrag der Amplitude berechnest, merkst du es natürlich nicht, doch die Phase ändert sich trotzdem ständig. Da der Goertzel nur eine vereinfachte FFT ist, macht sich eine Abtastverschiebung nicht im Betrag der Amplitude bemerkbar, sehr wohl jedoch in der Auswertung der Phase. Fazit: Bei der Abtastung ist eine kontinuierliche Verschiebung der Sampels verboten!
Wie kann ich mir das berechnete Spektum anschauen? Wären das Betrag und Phase nach jedem Schleifendurchlauf?
Mr. Qb schrieb: > Wie kann ich mir das berechnete Spektum anschauen? Wären das Betrag und > Phase nach jedem Schleifendurchlauf? Immer dann wenn ein kompletter Block berechnet ist.
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...
Wenn die Blocklänge N ist, so muss die Schleife von 0 bis N laufen. Der Wert x[N] muss Null sein.
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.
Angehängte Dateien:
-
Amplitude.png
4 KB -
Phase.png
4,6 KB -
real.png
2,9 KB -
imag.png
2,9 KB
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?
Mr. Qb schrieb: > aber meine Phase ist ein > Sägezahn. Einer ne Idee? Nun das sieht so aus, als ob du eine Phasenverschiebung von Block zu Block hast.
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?
Mr. Qb schrieb: > Kann ich den Filter auch als eine art Bandpass sehen? Das Filter realisiert exakt die FFT auf genau einer Frequemz mit der vorgegebenen Bandbreite.
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.
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.
Hi, ich habe auch noch eine Frage zum Goertzel. Ich möchte diesen zur Amplituden- und Phasenbestimmung nutzen, leider passt schon die Amplitude nicht wirklich und ich weiß nicht warum. Letztlich verwende ich den Code aus dem Thread:
1 | void SignalProcessing::testIt() |
2 | {
|
3 | //Signal generieren
|
4 | QList<double> sign; |
5 | int sampleRate=1000; |
6 | int periods=1; |
7 | int frequency=100; |
8 | int amp=1000; |
9 | double measuredAmplitude; |
10 | double measuredPhase; |
11 | int sampleCount=(sampleRate/frequency)*periods; |
12 | for (int i=0; i<sampleCount; i++) |
13 | {
|
14 | sign << amp*sin((2*M_PI*i)/(sampleRate/frequency)); |
15 | }
|
16 | goertzel(&test,frequency,sampleRate,&measuredAmplitude,&measuredPhase); |
17 | qDebug() << "Görzel:" << measuredAmplitude << measuredPhase; |
18 | }
|
19 | |
20 | void SignalProcessing::goertzel(QList<double> *inPutSignal, int frequency, int sampleRate, double *amplitude, double *phase) |
21 | {
|
22 | inPutSignal->append(0); |
23 | int N=inPutSignal->count(); |
24 | double D0=0; |
25 | double D1=0; |
26 | double D2=0; |
27 | double realPart; |
28 | double imagPart; |
29 | |
30 | double a1=2.0*cos(2.0*M_PI*frequency/sampleRate); |
31 | for(int i=0; i< N; i++) |
32 | {
|
33 | D0=inPutSignal->value(i)+a1*D1-D2; |
34 | D2=D1; |
35 | D1=D0; |
36 | }
|
37 | realPart=D1-(a1/2)*D2; |
38 | imagPart=D2*sin(2*M_PI*frequency); |
39 | *amplitude=sqrt(realPart*realPart+imagPart*imagPart); |
40 | *phase=atan(imagPart/realPart); |
41 | if (realPart<0) |
42 | {
|
43 | *phase=*phase-M_PI/2; |
44 | }
|
45 | else
|
46 | {
|
47 | *phase=*phase-M_PI/2; |
48 | }
|
49 | }
|
Für Amplitude kommt 0.000437114 und für die Phase kommt -1.5708 raus, was ja nicht stimmen kann. Hat jemand diesbezüglich eine Idee? Grüße MrPPI
Hallo Joe, vielen Dank für die schnelle Antwort. Ich generiere mir ja mit der ersten Funktion einen sinus (eine Periode). Mit
1 | inPutSignal->append(0); |
2 | int N=inPutSignal->count(); |
in der Goertzel-Funktion füge ich einen weiteren Punkt mit dem Wert 0 hinzu. Das sollte doch stimmen oder? Grüße MrPPI
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
So ich bin weitergekommen:
1 | void SignalProcessing::goertzel(QList<double> *inPutSignal, int frequency, int sampleRate, double *amplitude, double *phase) |
2 | {
|
3 | inPutSignal->append(0); |
4 | int N=inPutSignal->count(); |
5 | double D0=0; |
6 | double D1=0; |
7 | double D2=0; |
8 | double realPart; |
9 | double imagPart; |
10 | |
11 | double a1=2.0*cos(2.0*M_PI*frequency/sampleRate); |
12 | for(int i=0; i< N; i++) |
13 | {
|
14 | D0=inPutSignal->value(i)+a1*D1-D2;+a1*D1-D2; |
15 | qDebug() << inPutSignal->value(i); |
16 | D2=D1; |
17 | D1=D0; |
18 | }
|
19 | |
20 | realPart=D1-(a1/2)*D2; |
21 | imagPart=D2*sin(2*M_PI*frequency/sampleRate); |
22 | *amplitude=sqrt(realPart*realPart+imagPart*imagPart); |
23 | *phase=atan(imagPart/realPart); |
24 | if (realPart<0) |
25 | {
|
26 | *phase=*phase-M_PI/2; |
27 | }
|
28 | else
|
29 | {
|
30 | *phase=*phase-M_PI/2; |
31 | }
|
32 | }
|
Ich musste bei "imagePart=D2*sin(2*M_PI*frequency) das / sampleRate eintragen. Jetzt ist noch folgendes: Wenn ich mir mit dem Signalgenerator eine Periode des Signals genierer (z.B. mit Amplitude=10) kommt 50 raus. Bei 2 Perioden: 100; D.h. ich muss durch 10 teilen und die Anzahl der Perioden mit einrechnen, was ja aber eigentlich nicht sein kann. Irgendwo ist noch der Wurm drin. Grüße
Thread beobachten
Seitenaufteilung abschaltenBitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.