PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : Variometer Algorithmus



nichtgedacht
14.01.2019, 18:52
Hallo zusammen,

für ein Eigenbau Vario verwende ich zur Bestimmung der Steigrate das
Dual Filter Differenz Prinzip. Siehe auch
https://gregstanleyandassociates.com/whitepapers/FaultDiagnosis/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.



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

husi
16.01.2019, 23:12
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

nichtgedacht
17.01.2019, 08:34
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

nichtgedacht
21.01.2019, 19:42
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.com/whitepapers/FaultDiagnosis/Filtering/OtherFilters/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:


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