Netzmessung: Frequenz, Amplitude, Phase aus ADC-Samples

[Michael H. (overthere)]

Hallo,

folgendes Setup:
Ich möchte vom 50Hz Netz aus Frequenz, Amplitude & Phasenlage bestimmen. 
Das ganze mit möglichst wenig Rechenaufwand (=Laufzeit).
Cortex M4 habe ich zum Rechenen, also sollte einiges gehen.

Ich messe z.B. mit einer Abtastrate von 10kHz die Netzfrequenz und 
möchte (nach einer Einschwingzeit x) die oben genannten Werte erhalten.

Was ich mir bisjetzt überlegt hatte:
1. FFT
"Haste keine Ahnung mach 'ne FFT" und lies die entsprechenden Werte 
raus. Die Phasenlage würde ich aber so nicht bekommen.

Was mir suspekt ist: Viel zu viel Rechnerei.

2. Pony
Hier im Forum gefunden
Beitrag "Frequenz, Amplitude und Phase eines Sinussignals bestimmen"

Was mir suspekt ist: Abtastratensensibel

3. "PHASOR MEASUREMENT UNIT ALGORITHM USING DELAY DEVICES"
Irgendwie basierend auf Lisagou-Figuren, ohne das ich das jetzt richtig 
verstanden habe.
Demo
https://www.youtube.com/watch?v=ItjDd4kSHJc
Paper
http://www.eurasip.org/Proceedings/Eusipco/Eusipco2015/papers/1570093445.pdf

Sieht eigentlich ganz okay aus, jedoch immer noch Rechnerei.

4. Idee: Sinusfrequenz mit Bandpass filtern (z.B. 49Hz-51Hz) + 
Nulldurchgang messen

Welche Algorythmen ich schon ausgeschlossen habe:
-Nulldurchgangsmessung: Ungenau bei Netzmodulation


Meine Fragen an euch.
Habt ihr schon was implementiert, und wenn ja mit was habt ihr gute 
Erfahrungen? Wieviel Rechenzeit habt ihr einplanen müssen? Was sollte 
man auf jeden Fall bleiben lassen?

Danke für den Support,
Michael

[Hp M. (nachtmix)]

Michael H. schrieb:
> Irgendwie basierend auf Lisajous-Figuren, ohne das ich das jetzt richtig
> verstanden habe.

Ist aber das einfachste.
Ausser dem zu messenden Signal brauchst du eine genaue und stabile 
Referenzfrequenz, z.B. aus DCF77 oder GPS abgeleitet.

P.S.:
Fourier kannst du hier vergessen; der zeigt dir nur, ob das Signal 
verzerrt ist oder andere Frequenzen enthält.

[Mw E. (Firma: fritzler-avr.de) (fritzler)]

Frquenzmessung kannste ja schonmal so machen:
http://www.fritzler-avr.de/HP/netzfreq.php
oder mit einem internen 32Bit Zähler im µC.

Phasenlage vergleichste dann mit einem internen 50Hz Signal.

[Michael H. (overthere)]

Frizler, du machst das Quasi über einen gefilterten Nulldurchgang. Ich 
finde das mit den Transistoren eine clevere Lösung!

Ich habe die Hardware schon fertig (Netzmessung) und möchte eigentlich 
"nurnoch" einen Algo drüberprogrammieren.

[./. (Gast)]

> Phasenlage

ist ja wohl Kaese.
Woher soll die Bezugsphase kommen?
Von den Vogonen?

[Mw E. (Firma: fritzler-avr.de) (fritzler)]

./. schrieb:
>> Phasenlage
>
> ist ja wohl Kaese.
> Woher soll die Bezugsphase kommen?

Wie ich schon sagte, dazu braucht man "einfach" ein genaues internes 
50Hz Signal und vergleicht dazu.
Mir fällt aber nicht ein wozu man das bräuchte.

Michael H. schrieb:
> Ich habe die Hardware schon fertig (Netzmessung) und möchte eigentlich
> "nurnoch" einen Algo drüberprogrammieren.

Das Problem ist eben, dass die Netzfrequenz auf 10mHz genau gemessen 
werden muss.
Bei 10kS/s und ner FFT sehe ich nicht wie das klappen soll.
Eine FFT mit einem Binning <10mHz würde ich jetzt auch nicht als 
"möglichst wenig rechenaufwand" einstufen.
Meine Schaltung tastet das Signal quasi als 1Bit ADC mit 16MS/s ab.
Demnächst kommt dann noch das 1pps Signal aus nem GPS Empfänger dazu.

[Wolfgang H. (Firma: AknF) (wolfgang_horn)]

Hi, Michael,

das war meine erste Anwendung eines Switched-Capacitor-Filters von 
Reticon. Damals R5622, Nachfolger LMF100.
Beide Filter im IC gemeinsam als Bandpass geschaltet, RC-Tiefpässe davor 
und danach.
Sehr wenig Hardware für die Qualität der Filterflanken.

Ciao
Wolfgang Horn

[Lurchi (Gast)]

Die FFT bringt eher nicht die nötige Auflösung und ist auch nicht so 
wenig Rechenaufwand. Die Rechenzeit ist dabei eher weniger das Problem, 
sondern der Speicherbedarf kann auf einem µC zum Problem werden.

Die Phasor-methode sollte so oder so ähnlich wie im Artikel beschrieben 
schon ganz gut funktionieren. Bei den Sinus / cosinus werten kann man 
ggf. Tabellen nutzen, wenn die Rechenzeit knapp werden sollte. Der Witz 
bei der Methode ist, das nach dem ersten schritt mit einer sehr viel 
niedrigeren Zwischenfrequenz (nahe 0) gearbeitet wird. Entsprechend 
reichen da wenige Samples aus, und nur der erste Schritt ist wirklich 
zeitkritisch. Auch braucht man nur relativ wenig daten-speicher.

Die Filtermethode kann in Prinzip auch funktionieren, wenn der 
numerische Filter wirklich gut ist, aber dann wird auch der Aufwand für 
die Berechnung sehr hoch, höher als bei der Phasor methode.

[H.Joachim S. (crazyhorse)]

Phasenlage zu bestimmen verstehe ich trotzdem nicht. Macht doch 
keinerlei Sinn die Phasenlage zu einem beliebigen, von mir aus auch 
hochgenauen, aber dennoch asynchronen Referenzsignal zu messen? Solange 
die Signale nicht  phasenstarr sind, gibts da nichts zu messen. Da kann 
jeder beliebige Wert rauskommen, der dann auch noch ständig wandert, 
wenn die Signale nicht frequenzgleich sein (und das können sie nicht 
sein ohne Kopplung). Aber selbst bei Frequenzgleichheit sagt das doch 
überhaupt nichts aus?

[m.n. (Gast)]

H.Joachim S. schrieb:
> Solange
> die Signale nicht  phasenstarr sind, gibts da nichts zu messen.

Ich denke, das ist beim TO noch nicht angekommen.

[Reinhard M. (Gast)]

./. schrieb:
>> Phasenlage
>
> ist ja wohl Kaese.
> Woher soll die Bezugsphase kommen?

Dachte ich auch im ersten Augenblick.
Aber für die Energieversorgungsunternehmen kann es durchaus
vorteilhaft sein die Einspeisung in das Verbundnetz,
das dann als Referenz dient, frequenz- UND phasenrichtig durchzuführen.
Der Phasor, wie in dem obengenannten Paper beschrieben, ist dann
ein vernünftiger Messwert für die Regelung.

Das Paper weißt ja auch auf eine Anwendung bei einem japanischen EVU 
hin.
Für den Endanwender im trauten Heim ist die Phase natürlich beliebig 
sinnvoll ;-)

[Helmut L. (helmi1)]

Reinhard M. schrieb:
> Aber für die Energieversorgungsunternehmen kann es durchaus
> vorteilhaft sein die Einspeisung in das Verbundnetz,
> das dann als Referenz dient, frequenz- UND phasenrichtig durchzuführen.

Dann dient aber das 2. Netz als Bezugspunkt um die Phasengleichheit 
herzustellen. In den Anfaengen der Wechselstromtechnik wurde sowas mit 
Glimmlampen in Hell oder Dunkelschaltung syncronisiert.

[Reinhard M. (Gast)]

Um nochmal auf das Paper zu kommen:

Die Pseudo IQ Generierung finde ich clever und kann ich in
einer C++ Simulation erfolgreich nachvollziehen.

void IQ_gen (double x)  // 2.4 IQ-Signal Generation
{
    static double x_1 = 0.0;
    double I = x;
    double Q = a0 * x + a1 * x_1;
    x_1 = x;

    //fprintf(fp,"%f %f\n",I,Q);

    x1.clear();
    x1.push_back(I);
    x1.push_back(Q);
}

Der Referenzoszillator ist simpel

void REF_gen (double x_ref)
{
    double I = sin(x_ref);
    double Q = -cos(x_ref);

    x2.clear();
    x2.push_back(I);
    x2.push_back(Q);
}

Die Berechnung des Phasors liefert sinnvolle Resultate

void vector_op()  //2.3 Vector Operation
{
    double skalarprodukt = x1[0] * x2[0]  +  x1[1] * x2[1];
    double betrag_x1 = sqrt( x1[0]*x1[0] + x1[1]*x1[1] );
    double betrag_x2 = sqrt( x2[0]*x2[0] + x2[1]*x2[1] );
    double res = /*acos*/ ( skalarprodukt /( betrag_x1 * betrag_x2) );

    double determinante =  x1[1]*x2[0] - x1[0]*x2[1];
    double res2 = /*asin*/ ( determinante /( betrag_x1 * betrag_x2) );

    x3.clear();
    x3.push_back(res);
    x3.push_back(res2);

    //fprintf(fp,"%f \n",res);

}

Aber die eigentliche Frequenzmessung habe ich nicht verstanden.
Meine Interpretation, die NICHT funktioniert.

void F_meas(double s)
{
    static double SP[2] ={0.0,0.0};
    static double SQ[2] ={0.0,0.0};
    static double SR[2] ={0.0,0.0};
    static double SS[2] ={0.0,0.0};

    SP[0] = s;
    SQ[0] = SP[1];
    SR[0] = SQ[1];
    SS[0] = SR[1];

    double LS1 = SQ[0] * SR[1] - SQ[1] * SR[0];
    double LS2 = SP[0] * SS[1] - SP[1] * SS[0];

    SP[1] = SP[0];
    SQ[1] = SQ[0];
    SR[1] = SR[0];
    SS[1] = SS[0];

    double f = 
    (10000.0 / ( 2.0 * TWOPI)) *  (acos (  1.0/2.0 * LS2/LS1 -1) );

    
    cout <<  f << endl;
}

Hat sich schon jemand mit der Frequenzmessung im Detail beschäftigt ?

[Reinhard M. (Gast)]

Die Frequenzmessung klappt jetzt:

double Fmeas (double s_0, double fsample)   // 2.2 Frequency Measurement
{
  double K = 2.0;
  static double s_1 = 0.0;
  static double s_2 = 0.0;
  static double s_3 = 0.0;
  static double s_4 = 0.0;

  double LS1 = s_1 * s_2  -  s_0  * s_3;
  double LS2 = s_0 * s_3  -  s_1  * s_4;
  double f = (fsample/(TWOPI * K)) *acos(0.5*LS2/LS1);

  s_4 = s_3;
  s_3 = s_2;
  s_2 = s_1;
  s_1 = s_0;


  return f;


}

Zur Simulation lasse ich einen Tiefpassfilter mitlaufen,
um eine Phasenverschiebung zu generieren:

/***********************************************************************************************************************************\
              |--> DUT (filter,etc) --sin(nW + phi)-> PseudoIQ ---(sin(nW + phi),-cos(nW + phi))--->|-----------|
  sin(nW)---->|                                                                                     | Vector_op |----> phasor--> phi
              |-------------------------------------> PseudoIQ ---(sin(nW),-cos(nW))--------------->|-----------|
                                                     (Referenz)
\***********************************************************************************************************************************/
void analysePhasor()
{
    FILE *fp;
    fp = fopen( "sniff.txt", "w+");

    double fsample = 10000.0;
    double fsignal = 50.123456;
    double fcut = 50.0;
    double f=0;

    W = TWOPI * fsignal / fsample; // Ω
    rho = W - PI/2.0;

    a0 = tan(rho);
    a1 = 1.0 / cos(rho);


    // create a lpfilter
    Biquad * lp = new Biquad();
    lp->setBiquad( bq_type_lowpass, fcut/fsample , 0.707, 0);

    for (int n = 0; n <= 2000; n++)
    {
        // Signal
        double x = sin(  n * W );
        if(n <= 4)
      f = Fmeas(x,fsample);

       // Filter
        double xx = lp->process(x);


        // Pseudo IQ
        IQ_gen(xx);

        // Referenz
        REF_IQ_gen(x);

        // Vector Operation
        vector_op();
        double w_rad = acos(x3[0]);
        double w_deg = w_rad / PI * 180.0;

        // Print Result(s)
        printf("%f  %f    %f  %f\n",x,xx,w_deg,f);

        // gnuplot -p -e "plot 'sniff.txt' with lines"
         fprintf(fp,"%f %f\n",x,xx);            // Lissajou
        //fprintf(fp,"%f %f\n",x1[0],x1[1]);    // IQ
        //fprintf(fp,"%f\n",w_deg);             // phi
        //fprintf(fp,"%f\n",xx);                // Filter output

    }
    fclose(fp);
}

quick und dirty

 g++ main.cpp Biquad.cpp -o PMU_Simulation

[Michael H. (overthere)]

Hallo Reinhard,

ich bin beeindruckt.

Deinen Code muss ich gleich mal nachvollziehen.

-Michael

[Reinhard M. (Gast)]

Hallo Michael
über die F -Messung muss noch meditieren.
Es sieht so aus, als müsste der Beginn der Messung beim Nulldurchgang
starten.
Ein Rücksetzen der Zustandsvariablen bringt nichts.
deshalb

       if(n <= 4)
    f = Fmeas(x,fsample);

Naja, für eine erste Demo reicht es.
Ich möchte mich weiter mit dem Thema beschäftigen, falls ich die
Zeit dazu finde.

Schönes WE.

[Michael H. (overthere)]

Naja, wie Du schon geschrieben hast, die F-Messung will noch nicht; die 
misst nur das 4. mal richtig, wenn man es weiterlaufen lässt, dann will 
es nicht.
Ich werde mir mal das Paper nochmals genau durchlesen (müssen) und noch 
ein paar Grundlagen verstehen, denn soweit wie Du bin ich noch nicht.

[Michael H. (overthere)]

Update: Was mir bei der Frequenzmessung aufgefallen ist: Wir haben ja 
einen Delay von 4, also brauchen wir 5 variablen.

Der hat da ja eine Mittelwertbildung reingebracht, leider ist die Formel 
nicht so wirklich leicht verständlich.

Hat jemand das Lissajous Product verstanden, wie das funktioniert?

Was mir sonst noch bei der Frequenzmessung aufgefallen ist:
Wird das nicht bei jedem durchlauf zu Null gesetzt?
>static double s_1 = 0.0;
>static double s_2 = 0.0;
>static double s_3 = 0.0;
>static double s_4 = 0.0;

[Mw E. (Firma: fritzler-avr.de) (fritzler)]

Guck dir mal unbedingt an was "static" in C bedeutet.

[Michael H. (overthere)]

Static: Ja, habs nochmal nachgelesen. Passt scho. Hatte ich falsch in 
Erinnerung...
Variablenanzahl: Auch klar, wir haben die ganze Zeit schon 5.

Jedoch sagt, mir meine Nase, das in diesem L.-Produkt der Hund vergraben 
liegt. Im Paper ist das nicht gut erklärt.

[Horst (Gast)]

Michael H. schrieb:
> 1. FFT
> "Haste keine Ahnung mach 'ne FFT" und lies die entsprechenden Werte
> raus. Die Phasenlage würde ich aber so nicht bekommen.

Warum das denn nicht? Natürlich liefert die FFT auch die Phasenlage.
Und wenn eine hohe Auflösung bei nur wenigen interessanten Bins zu 
aufwändig ist, dann schaue man sich den Görtzel-Algo an.

[Michael H. (overthere)]

So, ich hab die Frequenzmessung am laufen.

Ich habe die alten s-Werte in den Kommentaren stehen, die neuen, 
richtigen im Code.

double s_1 = 0.0;
double s_2 = 0.0;
double s_3 = 0.0;
double s_4 = 0.0;
double s_5 = 0.0;


double Fmeas (double s_0, double fsample)   // 2.2 Frequency Measurement
{
  double K = 2.0;

  double LS1; //= s_1 * s_2  -  s_0  * s_3;
  double LS2; // = s_0 * s_3  -  s_1  * s_4;
  
  LS1=s_1*s_4-s_2*s_3;
  LS2=s_0*s_5-s_1*s_4;
  
  double f = (fsample/(TWOPI * K)) *acos(LS2/(2*LS1));

  fprintf(fp,"%10f %10.8f %10.8f %10f %10f %10f %10f %10f\n",f, LS1, LS2, s_0,s_1,s_2,s_3,s_4);
  
  s_5 = s_4;
  s_4 = s_3;
  s_3 = s_2;
  s_2 = s_1;
  s_1 = s_0;

  return f;


}

Was jetzt noch ein "lustiges" Phänomen ist:
Selbst eine kleine Störung 10^-6 Peak-RMS liefert große Abweichungen.

double x = sin(  n * W ) + (rand()%100)*0.00000001;
double xref = sin(  n * W_ref);
        
f = Fmeas(x,fsample);

Erzeugt mir folgenden Sniff:

      -nan 0.00000000 0.00000000   0.000001   0.000000   0.000000   0.000000   0.000000
      -nan 0.00000000 0.00000000   0.031489   0.000001   0.000000   0.000000   0.000000
      -nan 0.00000000 0.00000000   0.062946   0.031489   0.000001   0.000000   0.000000
1250.000000 -0.00000003 0.00000000   0.094340   0.062946   0.031489   0.000001   0.000000
1249.984281 -0.00198204 -0.00000008   0.125642   0.094340   0.062946   0.031489   0.000001
 49.747421 -0.00198199 -0.00395623   0.156818   0.125642   0.094340   0.062946   0.031489
 50.039151 -0.00198201 -0.00395618   0.187839   0.156818   0.125642   0.094340   0.062946
 51.001307 -0.00198212 -0.00395611   0.218674   0.187839   0.156818   0.125642   0.094340
 48.596202 -0.00198192 -0.00395646   0.249291   0.218674   0.187839   0.156818   0.125642
 51.740191 -0.00198224 -0.00395610   0.279662   0.249291   0.218674   0.187839   0.156818
 50.246238 -0.00198210 -0.00395630   0.309755   0.279662   0.249291   0.218674   0.187839
 48.372951 -0.00198181 -0.00395630   0.339541   0.309755   0.279662   0.249291   0.218674
 51.109082 -0.00198205 -0.00395592   0.368991   0.339541   0.309755   0.279662   0.249291
 49.941417 -0.00198201 -0.00395621   0.398074   0.368991   0.339541   0.309755   0.279662
 49.441582 -0.00198202 -0.00395639   0.426762   0.398074   0.368991   0.339541   0.309755
 51.774282 -0.00198229 -0.00395619   0.455027   0.426762   0.398074   0.368991   0.339541
 49.335792 -0.00198204 -0.00395646   0.482841   0.455027   0.426762   0.398074   0.368991
 50.181511 -0.00198200 -0.00395611   0.510176   0.482841   0.455027   0.426762   0.398074
 49.873847 -0.00198192 -0.00395605   0.537006   0.510176   0.482841   0.455027   0.426762
 49.373295 -0.00198190 -0.00395618   0.563302   0.537006   0.510176   0.482841   0.455027
 50.184781 -0.00198215 -0.00395642   0.589040   0.563302   0.537006   0.510176   0.482841
 53.642886 -0.00198243 -0.00395586   0.614194   0.589040   0.563302   0.537006   0.510176

Folglich: Filtern ist Pflicht!

Reinhard: Wie würdest Du das machen?

[Michael H. (overthere)]

So, ich habe mal einen universalen Filter implementiert, mit etwas mehr 
Rechenaufwand (leider natürlich...)

#define NUM_OF_AVG 50
#define NUM_OF_MEM (5+NUM_OF_AVG)

double sample[NUM_OF_MEM];
double Fmeas_uniavg (double s_0, double fsample)   // 2.2 Frequency Measurement
{
  double K = 2.0;

  double LS1=0;
  double LS2=0;
  sample[0]=s_0;
  int i;
  for(i=0;i<NUM_OF_AVG;i++){
    LS1=LS1+(sample[0+i]*sample[3+i]-sample[1+i]*sample[2+i]);
    LS2=LS2+(sample[0+i]*sample[5+i]-sample[1+i]*sample[4+i]);
  }
  
  double f = (fsample/(TWOPI * K)) *acos(LS2/(2*LS1));

  fprintf(fp,"%10f %10.8f %10.8f %10.8f %10f %10f %10f %10f %10f\n",f, (LS2/(2*LS1)), LS1, LS2, sample[0],sample[1],sample[2],sample[3],sample[4]);
  
  for(i=(NUM_OF_MEM-1);i>0;i--){
    sample[i]=sample[i-1];
  }
  
  return f;


}

Einen Kommentar noch zur Genauigkeit:
Je kleiner man den Änderungswinkel macht, desto mehr wirkt sich das 
Rauschen aus. Daher ist es empfehlenswert die Samplingfrequenz zu 
senken, und eventuell mit dem MCU Hardware-Averaging betreiben.

Wenn ich jetzt den Filter auf 50 Durchschnittswerte und die 
Samplingfrequenz auf 1000Hz Einstelle komme ich auf 50.1234xx was für 
OTTO-ADVANCED-Programmer reichen sollte.

Zum Schluss muss ich noch was loswerden: Der Typ ist ein Genie. Ich 
staune immer wieder was man alles mit Sinus/Cos-Sätzen machen kann.

[Mw E. (Firma: fritzler-avr.de) (fritzler)]

Ist jetzt nurnoch die Frage ob das Ergebnis stimmt.
Kannste ja hiermit abgleichn: http://www.netzfrequenzmessung.de/
Aber nicht daurnd laufen haben, nach ner Weile wird man für nen Tag 
gesperrt wegen zu viel Traffic.
Server zu schwach...

[Michael H. (overthere)]

Mw E. schrieb:
> Ist jetzt nurnoch die Frage ob das Ergebnis stimmt.

Ich erzeuge im Code meine eigene Referenzfrequenz, Das Skript schiebt 
nur Zahlen hin und her. Und daher weiß ich, dass es stimmt. Weil das 
rauskommt was rauskommen soll.

[Reinhard M. (Gast)]

Hallo Michael
entschuldige das ich mich erst jetzt wieder melde,
aber am Wochenende bin ich meistens nicht online;-)

Du hast vollkommen recht: Filtern ist Pflicht.
In dem Paper wurde dem Averaging auch ein Absatz gewidmet,
den hab ich bislang ignoriert, weil ich ihn nicht so wichtig hielt.
Dem ist aber anscheinend nicht so.

Deine Updates werde ich mir dann gleich mal anschauen.

Michael H. schrieb:
> Je kleiner man den Änderungswinkel macht, desto mehr wirkt sich das
> Rauschen aus. Daher ist es empfehlenswert die Samplingfrequenz zu
> senken, und eventuell mit dem MCU Hardware-Averaging betreiben.
>
> Wenn ich jetzt den Filter auf 50 Durchschnittswerte und die
> Samplingfrequenz auf 1000Hz Einstelle komme ich auf 50.1234xx was für
> OTTO-ADVANCED-Programmer reichen sollte.

Das ist schon ein sehr interessantes Phänomen.

[Michael H. (overthere)]

Ich habe an der Schaltung etwas weiter gearbeitet und auch mal ein 
GitHub Repo erstellt.

https://github.com/greenscreenflicker/grid-meas