• Liebe RC-Network Benutzer: Bitte beachtet, dass im August 2020 alle Passwörter zurückgesetzt wurden. Mehr dazu in den News...

Variometer Algorithmus

Hallo zusammen,

für ein Eigenbau Vario verwende ich zur Bestimmung der Steigrate das
Dual Filter Differenz Prinzip. Siehe auch
https://gregstanleyandassociates.co...sis/Filtering/MACD-approach/macd-approach.htm

Dabei wird die Eigenschaft eines Tiefpassfilters nützlich, sich der Steigung des Eingangssignals anzunähern und dieser
dann zeitlich um den Betrag Tau verzögert zu folgen. Die Differenz der Ausgänge zweier paralleler Tiefpassfilter ist dabei
der Steigung Proportional. Siehe Bild 1. In dem gelben rechtwinklige Dreieck ist die Steigung gebildet aus der Differenz der
Ausgangssignale geteilt durch die Differenz der Zeitkonstanten und entspricht genau der Steigung der blauen Linie
wenn der Vorgang einige Tau lang andauert.

Bei passender Wahl der Zeitkonstanten erhält man so die Erste Ableitung der Kurve (Steigung) bei gleichzeitiger Unterdrückung
des Rauschens der Sensordaten. Während aber das gewonnene Signal für die Höhe bereits durch den Tiefpassfilter mit der
kleineren Zeitkonstanten gut geglättet erscheint, ist das Signal für das Steigen nicht so gut geglättet. Das liegt daran, dass
die Ableitung der Steigung eines Verlaufs ja mit einem wie auch immer gearteten Hochpassfilter erfolgen muss. Und ein
Hochpassfilter verstärkt ja prinzipiell das Rauschen.

Es liegt nun also nahe das Variometersignal selbst durch einen weiteren Tiefpass zu schicken.
Dabei tritt aber eine weitere Verzögerung auf.

Was man bei einem Vario haben will ist ja eine möglichst spontane Reaktion ohne zappelig zu sein.

Das Rauschen des MS5611 Sensors ist vorhersehbar und bleibt dauerhaft unter einer bestimmten Schwelle.
Bei schnellen Anstiegen des Signals wird die Kurve durch das Rauschen relativ weniger stark verzerrt.

In Anlehnung an https://onlinelibrary.wiley.com/doi/abs/10.1002/aic.690260120
habe ich den Fillterkoeffizienten dynamisch gemacht um diesen Umstand zu berücksichtigen.
Der Filter hat dadurch eine höhere Reaktionsgeschwindigkeit ohne die Glättung schlechter
werden zu lassen. Bild 2 und 3 zeigen den Unterschied. Das originale Sensorsignal
ist mit einem künstlich erzeugtem Höhenverlauf beaufschlagt, Siehe Code Block.

Code:
    if (ms5611.data_ready) {    // flag is interrupt trigggered and reset by ms5611.getPressure  

        long realPressure = ms5611.getPressure (true);
        double relativeAltitude = ms5611.getAltitude (realPressure, referencePressure);

        j++;
        if ( j <= 200) {
            relativeAltitude += j * 0.0125;  //2.5 m / 2.5 s
        }
        else {
            if ( j <= 280 ) {
                relativeAltitude += 2.5;     //2.5 m
            } else {
                if ( j <= 300 ) {
                    relativeAltitude += 2.5 - ( j - 280 ) * 0.1250; //2.5 m / 0.25 s
                }  
            }
        }    
 
        if ( j == 360) {
            j = 0;
        }
        
        r_altitude0 = r_altitude0 - 0.076923  * (r_altitude0 - relativeAltitude); // dt 12500 us T = 0.15
        r_altitude = r_altitude -  0.04 * (r_altitude - relativeAltitude);         //             T = 0.3
        climb0 = (r_altitude0 - r_altitude) * 6.66667;   // dt 12500 us
        
        dyn_alfa = abs( climb - climb0 / 1 );
        
        if ( dyn_alfa >= 1 ) {

            dyn_alfa = 1;
           
        }
                
        climb = climb - dyn_alfa * ( climb - climb0 );

        ...

Gruß
Dieter
 

Anhänge

husi

User
Hallo Dieter,

vielen Danke für die Einblicke in dein Projekt.
Das Thema finde ich sehr spannend, auch weil ich mich in mittlerer Zukunft damit beschäftigen möchte.

Kannst du mir bitte ein paar Fragen zu deinem Quellcode erläutern?
Du hast die Variable "j", die ja von 0 bis 360 hochgezählt wird und in Stufen (0 bis 199) -> (200 bis 279) -> (280 bis 360) den Additions-Therm (relative Höhe) verändert. Wofür steht die Schleife von 0 bis 360 und in welcher Zeit wird sie durchlaufen? Und warum bildest du diese flache "Rampe" mit einem kurzen geraden "Rücken", der am Ende steil bis auf 0 wieder abfällt?

Woher kommen die Zahlen 0,076923 bzw. 0,04 die anscheinend deine beiden Tiefpass-Filter sind?
Warum rechnest du die Differenz (climb - climb0) zwei mal aus und pufferst die Differenz nicht in einer eigenen Variablen? Ist der Microcontroller schnell genug und du hast gleichzeitig Platznot im Arbeitsspeicher?

Es würde mich sehr freuen, wenn du dir die Zeit für weitere Erklärungen nehmen würdest.


Vielen Dank
Mirko
 
Hallo Mirko,

die Abtastzeit ist 12,5 ms. Siehe die modifizierte Arduino Lib https://github.com/nichtgedacht/Arduino-MS5611-Interrupt.
Die aufsteigende Rampe ist 2,5 m in 2,5 s und die abfallende Rampe 2,5 m in 0,25 s.

Die Glättung habe ich experimentell ermittelt und bin auf Tau1 = 0,3 s und Tau2 = 0,15 s raus gekommen
Die Tiefpassfilter sind die zeitdiskrete Form eines einfachen RC Filters.

Formeln:

Tau = R * C
Alfa = delta_t / ( Tau + delta_t )

Steigung = dx / ( Tau1 - Tau2 )
Steigung = dx / ( (delta_t / Alfa1) - ( delta_t / Alfa2 ) ) // dx ist die Differenz der Tiefpassausgänge

Gruß
Dieter
 
Hallo zusammen,

ich habe noch weiter experimentiert.

Bei einem "double exponential filter" kann der Ausgang
des Filters dem Rampenverlauf (Ansteigen) des Eingangssignals
prinzipiell verzögerungsfrei folgen.
https://gregstanleyandassociates.co...otherfilters.htm#Double-Exponential-Smoothing

Durch Steuerung der Glättungszeitkonstanten
abhängig von der Steigung (siehe erstes Posting) kann man auch hier
die Glättung noch optimieren.

Die Idee:
die zwei parallelen Tiefpassfilter zur Gewinnung der Ableitung
(Vario) auf dem Ausgang von o.g. Filter arbeiten lassen und
wegen der quasi verzögerungsfreien Vorglättung des Rauschens
kleinere Zeitkonstanten wählen können.

Da es bei dem gesamten Vorhaben ( Vario ) nur um die Optimierung
im Zeitbereich geht kommt man auch durch Experimentieren im Zeitbereich
am schnellsten zum Optimum.

Man muss mit dem ersten Filter bei der Glättung einen Kompromiss
mit einem Überschwingen im Zeitbereich machen. Mit größerer Zeitkonstante
kann der Ausgang schnellen Änderungen der Steigung nicht folgen und der
Graph geht "schleudernd" auf die Mitte vom Rauschen zurück. Mit der
Zeitkonstanten der Glättung der injizierten Trend Variable kann man die
Frequenz dieses Überschwingens beeinflussen. Hier wurden beide
Zeitkonstanten relativ groß gewählt. Aus den Graphen ist eine
spontanere Reaktion des Varios zu erkennen. Weitere Optimierung ist
sicher möglich.

Die Linien im Diagramm:

Blau: Original Daten
Rot: Der "double exponential filter" Ausgang mit Steuerung
der Glättung abhängig von der Steigung.
Gelb: Der auf Rot aufsetzende 1. "exponential filter"
Grün: Der auf Rot aufsetzende 2. "exponential filter"
Violett: Die Höhendifferenz zwischen grün und gelb
Grau: Die mit einem "exponential filter" nachgefiltererte
o.g. Höhendifferenz mit Steuerung der Glättung abhängig
von der Steigung.

Code:
Code:
    if (ms5611.data_ready) {    // flag is interrupt trigggered and reset by ms5611.getPressure  

        long realPressure = ms5611.getPressure (true);
        double relativeAltitude = ms5611.getAltitude (realPressure, referencePressure);

        //relativeAltitude = 0;

        j++;
        if ( j <= 200) {
            relativeAltitude += j * 0.0125;  //2.5 m / 2.5 s
        }
        else {
            if ( j <= 280 ) {
                relativeAltitude += 2.5;     //2.5 m
            } else {
                if ( j <= 300 ) {
                    relativeAltitude += 2.5 - ( j - 280 ) * 0.1250; //2.5 m / 0.25 s
                }  
            }
        }    
 
        if ( j == 360) {
            j = 0;
        }


        dyn_alfa = abs( ( r_altitude0 - relativeAltitude ) / 2 );
        
        if ( dyn_alfa >= 1 ) {

            dyn_alfa = 1;
           
        }

        r_altitude0 = relativeAltitude + ( 1 - dyn_alfa ) * ( r_altitude0 + slope0 - relativeAltitude);

        slope0 = slope0 - 0.01 * ( slope0 - r_altitude0 + r_altitude0_prev );
 
        r_altitude0_prev = r_altitude0;

        r_altitude1 = r_altitude1 -  0.14285714 * (r_altitude1 - r_altitude0);
        r_altitude  = r_altitude - 0.076923  * (r_altitude - r_altitude0);
        
        climb0 = (r_altitude1 - r_altitude) * 12.6; 


        dyn_alfa = abs( (climb - climb0) / 2 );
        
        if ( dyn_alfa >= 1 ) {

            dyn_alfa = 1;
           
        }
                
        climb = climb - dyn_alfa * ( climb - climb0 );

        ...

Wenn jemand vergleichbare oder bessere Lösungen kennt bitte hier
zeigen.

Gruß
Dieter
 

Anhänge

mha1

User
Hallo Dieter,

habe das mit einem 8Mhz Arduino Mini nachgebaut. Ist natürlich besser als meine Versuche, d.h. ich passe mit besseren Lösungen.

Eines ist mir in der ms5611.cpp aufgefallen, weil sich der Compiler ueber eine Division durch Null beschwert:
TEMP2 = (dT * dT) / (2 << 30);

Die Konstante 2 wird als integer interpretiert. 30 mal nach links geschoben generiert das einen Überlauf. Ein simpler typecast behebt das:
TEMP2 = (dT * dT) / ((uint_32)2 << 30);

Lass uns wissen, wenn Du noch weitere Optimierungen Deiner Filter gefunden hast.

VG Michael
 
Servus Dieter,

ich hab mir auch ein Vario, allerdings zum Gleitschirmfliegen gebaut.
Zur Berechnung des Steigwertes hab ich einen PT2 Tiefpassfilter einmal abgeleitet.
Das hat den Vorteil, dass der ganze Algorithmus ohne eine Differenzbildung, sondern nur mit Integratoren aufgebaut werden kann.

Zum Code:
Der Filter kann mit den Konstanten d (Dämpfung) und T (Zeitkonstante) eingestellt werden.
Für mich haben sich beim Gleitschirmfliegen folgende Werte bewährt:

Code:
float T  	= 1;
float d 	= 0.9;
Dabei wird schön geglättet und auch bei recht turbulenten Bedingungen nicht jedes winzigste Heberchen angezeigt (passiert eh über´s Popometer).



Code:
// ***** PT2D Filter *****
	u = df_data->height;

	// Init Filter
	if (flagfirst)
	{
		// Init Filter
		yi1 = u;
		ui1 = u * ts;
		yi2 = 2 * ui1 - 2 * d * T * yi1;

		flagfirst = 0;
	}

        // Filter
	sub = ui1;
	ui1 += u * ts - sub;
	yi2 -= sub;

	y = (ui1 - 2 * d * T * yi1 - yi2) / (T * T);

	yi1 += y * ts;
	yi2 += yi1 * ts;

Noch eine kurze Erklärung zum folgenden Absatz:
Code:
	sub = ui1;
	ui1 += u * ts - sub;
	yi2 -= sub;

ui1 und yi2 werden bei jedem Durchlaufen immer größer, wodurch die Genauigkeit des Ergebnisses deutlich kleiner wird, da irgendwann auch ein Float nicht mehr reicht, klein genug aufzulösen.

Durch die Subtraktion von ui1 bleibt ui1 immer nahe 0. Durch die Subtraktion von ui1 von yi2 wird das Gesamtergebnis nicht verfälscht und auch yi2 bleibt nahe einem festen Wert.


Wäre cool, wenn du die genannten Filter mal alle implementierst und miteinander vergleichst ;)

Schöne Grüße,
Jakob
 

hul

User
das Dual Filter Differenz Prinzip hab ich bei der Arbeit schon verwendet (in Kraftwerksleittechnik).

Das mit dem einmal abgeleiteten PT2 ist mir neu. Jakob, kannst Du das bitte ein wenig erläutern? Ich scheitere bei dem Code.

Besten Dank im Voraus, Hans
 

mha1

User
Hallo,

ich habe kurz eine Simulation gebastelt um vergleichen zu können.

Das simulierte Eingangsignal ist ein ruhender Sensor, d.h. er liefert konstat 0m relative Höhe. Darauf ist ein zufälliges Rauschen mit +/-0.3m aufgeschaltet und darauf dann Dieters Rampe mit Steigen auf 2.5m in 2.5s, Halten der Höhe 2.5m für 1s, Sinken auf 0m in 0.25s und Bleiben auf 0m für 0.75s. Die Rampe dauert also 4.5s. Dann beide Filter (Dieters aus Post #4 und Jakobs aus Post #6) implementiert. Dieters bekommt die Daten alle 12.5ms abgetastet, weil die Zeitkonstanten wie er schreibt für 12.5ms bestimmt wurden. Jakobs Filter wird mangels besseres Wissens alle 20ms gerechnet. Das Eingangssignal und die beiden Filteroutputs werden gelogged und araus ein Excel Chart gemalt.

Wenn ich alles richtig gemacht habe, steht es 1:0 für Dieter. Vor allem das nahezu latenzfreie Folgen mit guter Rauschunterdrückung ist echt gut. Jakobs Filter entrauscht auch gut, genehmigt sich aber schon deutlich mehr Zeit dem Eingangssigmal zu folgen. Ich denke Jakobs Filter ist sicher für das deutlich trägere manntragende System entworfen worden.

Ein wenig Hilfe brauche ich auch noch. Woher kommt der Faktor 12.6 in Dieters Filter. Die Taus konnte ich mit 0.081s und 0.1561s und damit die (1-a) Faktoren gut nachrechnen. Aber die 12.6? Sieht für mich wie eine Art Skalierungsfaktor aus. Und in Jakobs Filter habe ich die Variable ts als sample Zeit interpretiert. Habe ich mit 20ms angenommen.

Dieters Filter habe ich auch in einem Arduino Mini Pro 8Mhz laufen, musste allerdings die ms5611 Bibliothek anpassen. Die Conversion mit einer Interruptroutine zu steuern ist eine gute Idee, aber wenn man Softserial verwenden möchten nicht machbar. Softserial nutzt auch wie die wire library pin change interrupts. Und das beisst sich, deshalb umgebaut. Damit sind aber die 12.5ms Zykluszeit nicht mehr ganz erreichbar, d.h. ich habe auf 50Hz umgebaut und dafür die Taus angepasst.

Viele Grüße

Michael

vergleich.png
 

Ay3.14

User
Interessantes Thema!

Der MS5611 ist, im Vergleich zum BMP180 oder BMP280, ein recht brauchbarer Luftdruck-Sensor (im Hobbybereich).

Schade ist allerdings, das (immer noch) eine Nachbearbeitung der Sensordaten erforderlich ist. Aber vielleicht ändert sich das auch irgendwann mal zum Besseren.

Danke für die Informationen.
 
Schade ist allerdings, das (immer noch) eine Nachbearbeitung der Sensordaten erforderlich ist. Aber vielleicht ändert sich das auch irgendwann mal zum Besseren.
Moin,

letztlich bestimmen die Anforderungen, wie man die rauschbehafteten Einzeldaten verarbeiten will.
Präzisionsbarometer oder schnelles Vario sind die Gegensätze.
Solange das Rauschen > Null ist will man als Anwender doch selber die Kontrolle haben.

Gruß
Dieter
 
mit einem analogen Drucksensor wie dem MPXA6115A sind eine schnelles und präzises Vario keine Gegensätze.
Moin

Kannst Du das mal näher erläutern?

Rauscht dieser Sensor gar nicht?
Welche Vorteile bezüglich Rauschen soll hier eine analoge Filterung vor der AD-Wandlung haben gegenüber erst AD zu wandeln und dann digital zu filtern?
Was ist im ersteren Fall mit dem Rauschen, das die abschließende AD Wandlung einbringt?

Wenn dann also letztlich doch Rauschen da ist, macht es doch für ein Barometer Sinn, das Rauschen einfach durch Mittelwertbildung zu entfernen. Die Werte werden ja höchstens alle paar Sekunden gebraucht.
Dafür kann man aber große Genauigkeit erreichen, wenn man auch noch die Temperatur stabilisiert.

Für ein schnelles und rauscharmes Variometer braucht man bessere Filter, was ja hier mein Thema ist.

Gruß
Dieter
 
Hallo Dieter,

bei dem analogen Drucksensor besitzt das Sensorsignal ein hohe Stabilität. Das Sensorsignal (0.2 - 4.8 Volt) wird über eine Schaltung mit Operationsverstärkern verstärkt, gefiltert, differenziert und skaliert. Die Berechnung der Steiggeschwindigkeit erfolgt also analog. Dieses Variosignal wird dann A/D gewandelt. Die Ansprechzeitkonstant verglichen zu den Varios mit MS5611 ist deutlich kleiner. Solche analogen Varios wie die alten wsTech sind natürlich Schaltungstechnisch wessentlich aufwendiger und daher auch entsprechend teuer. Dafür sind solche Varios sehr schnell und sie besitzen eine hohe Emfindlichkeit im Bereich von wenigen cm/sec. Damit kann auch kleinste Thermik noch schnell und präzise erfasst werden.

Ich hab mir vor Jahren ein Vario mit dem analogen Drucksensor MPX4115A gebaut. Das A/D gewandelte Variosignal sende ich über einen PIC 12F683an die serielle Schnittstelle eines FrSky D8 Empfängers. Mein Vario Dashboard wandelt das dann in Ton und Displayanzeige um.

Link zu meinem Vario, die erste Version: https://www.rcgroups.com/forums/showpost.php?p=17959296&postcount=93

Gruss
Micha
 

Ay3.14

User
"Rauschen" und Signalschwanken ist immer so eine Sache, aber wie ist es mit der "Sprungantwort" vom Sensor?

Wenn man zwei (oder mehr) Sensoren parallel (ort- und zeitgleich) auswerten würde, wäre dann das Ergebnis-Signal besser?
 
...
Die Ansprechzeitkonstant verglichen zu den Varios mit MS5611 ist deutlich kleiner.
...
Moin,

wenn Du hier von Deinem Analogfilter ein Schaltbild zeigen würdest, könnte man das in LTSPICE nachbilden und die Ergebnisse direkt vergleichen.
Ansonsten halte ich solche Pauschalaussagen für Unsinn.

@Ay3.14:

Auf den Ersten Blick könnte man meinen damit ein besseres Gesamtergebnis zu erreichen.
Z.B. Immer nur den Wert benutzen, der näher am aktuellen gefilterten Wert liegt.

Ich würde intuitiv aber lieber alternierend auslesen.

Gruß
Dieter
 

husi

User
Hallo Micha,

da klingt sehr gut und auf jeden Fall interessant. Wäre es für dich möglich, das Vario mit dem analogen Sensor und ein "normales Vario" (z.B. den GPS-Logger 3) gleichzeitig 2 oder 3 Etagen eine Treppe hoch und wieder runter zu tragen und uns beide Logs als Grafik oder Liste zu zeigen?

Das würde mich sehr freuen!

Viele Grüße
Mirko
 
Hallo Mirko,

ich besitze neben meinem DIY-Vario nur 5 alte wsTech Varios. Zudem kann ich die Daten von meinem DIY-Vario nicht loggen. Einen Vergleich kann ich daher nicht machen. Der Herr Schreiner (wsTech) hat das mal gemacht. Siehe http://wstech.de/ unter Variometerkunde, das Bild Treppenlauf mit Variometern.

Ich bin im Herbst oft mit meinem Thermik XL bei sehr schwacher Thermik geflogen. Aufwinde mit einem Steigen im Bereich von nur 0.3 m/sec werden vom eingebauten Voice Vario präzise, stabil und schnell erfasst. Die Varios mit Drucksensoren mit digitaler Sensorschnittstelle kommen an meine wsTech Varios nicht heran. So ist mein Eindruck im Vergleich zu normalen Varios einiger Vereinskollegen.

@ Dieter

Im Link hat es die Schaltung als PDF von meinem Vario, die endgültige Schaltung hab ich leider nicht mehr Dokumentiert. Da ich von meinem Vater vor einigen Jahren noch 3 wsTech Varios übernommen habe, hab ich nicht mehr am DIY-Vario weitergemacht. Die Mischung aus Analog-und Digitaltechnik mit SW ist zudem auch sehr Zeitaufwendig.

https://www.rcgroups.com/forums/showpost.php?p=26026220&postcount=1385

Gruss
Micha
 
Oben Unten