Effizienter Running Median

[Christoph M. (mchris)]

Ein Medianfilter ist gut geeignet, um Ausreißer wegzufiltern.

Ich bin auf der Suche nach einer Effizienten Routine in C++.

Es gibt scheinbar haufenweise Implementierungen, aber welche Art ist die 
effizienteste?

Im Moment scheint mir die Verwendung einer "linked List" am besten:
https://community.infineon.com/gfawx74859/attachments/gfawx74859/PSoC6MCU/8068/2/MedianFilter_v0_0_A.pdf

Gibt es schnellere Möglichkeiten?

[Niklas G. (erlkoenig)]

Da scheint eine lineare Suche beim Hinzufügen nötig zu sein? Da würde 
ich doch lieber eine sortierte Baumstruktur nehmen, std::set oder so.

[Katrin I. (Gast)]

Kommt auf die Menge der zu erfassenden Daten und den Bedarf an Rückschau 
an. Zu beiden hat sich der TE noch nicht eingelassen.

[Christoph M. (mchris)]

Die Ausreißer sind eher selten und der Filter soll eine Länge zwischen 
3,5,7,9 Samples haben.

Für den Rechenzeitbedarf würde ich kleiner gleich 20us auf einem 
Atmega328 für die Länge 3 und int16_t erwarten.

[Norbert (der_norbert)]

Christoph M. schrieb:
> Für den Rechenzeitbedarf würde ich kleiner gleich 20us auf einem
> Atmega328 für die Länge 3 und int16_t erwarten

Das wären mehr als reichliche 320 Taktzyklen bei 16MHz.
Da kann man's ja fast in COBOL schreiben.

[Michael B. (laberkopp)]

Christoph M. schrieb:
> Gibt es schnellere Möglichkeiten?

It depends...


Bei kleinen Filterlängen, 3,5,7, ist binäre Suche in einem Array dessen 
Elemente verschoben werden oft schneller. Bei grossen Filterlängen ein 
B-Baum.

[Ob S. (Firma: 1984now) (observer)]

Michael B. schrieb:

> It depends...

Absolut.

> Bei kleinen Filterlängen, 3,5,7, ist binäre Suche in einem Array dessen
> Elemente verschoben werden oft schneller. Bei grossen Filterlängen ein
> B-Baum.

Die genaue Grenze hängt von der Architektur ab, vom Compiler und vor 
allem auch vom Signal selber. Es ist also sehr schwierig bis nahezu 
unmöglich, hier eine genaue genaue Grenze für die Entscheidung Array vs. 
B-Tree anzugeben. Da hilft nur: beides Implementieren und dann 
Ausprobieren, was im konketen Fall besser performt.

Und nicht mal dann ist man auf der sicheren Seite: die 
Signaleigenschaften können sich schließlich auch ändern...

[Katrin I. (Gast)]

Michael B. schrieb:
> Bei kleinen Filterlängen, 3,5,7,

Was heißt denn hier Filterlänge? Zahl der zu untersuchenden letzten 
Samples? Wieviele sollen behalten werden? Und wie wird der Median 
bestimmt? Warum nicht gleich addieren und den kleinsten / größten 
mitnehmen und jeweils entscheiden? Das wäre das Einfachste, wenn immer 
nur ein Ausreisser gefiltert werden muss.

Und dann braucht es noch eine Definition für Ausreisser, weil das ja 
auch ein sich ändernder Messwert sein kann.

[Christoph M. (mchris)]

>Was heißt denn hier Filterlänge?

https://de.mathworks.com/help/matlab/ref/movmedian.html

[Norbert (der_norbert)]

Kati schrieb:
> Und wie wird der Median bestimmt?

Das ist definiert und bereits millionenfach beschrieben.

> Warum nicht gleich addieren und den kleinsten / größten
> mitnehmen und jeweils entscheiden?

Das ist mit beinahe chirurgischer Präzision das exakte Gegenteil.

[Rbx (rcx)]

Christoph M. schrieb:
> Gibt es schnellere Möglichkeiten?

Möglicherweise C oder Asm. Wäre ich aber nicht so sicher, die 
verfügbaren Bibs sind schon sehr gut.
Den Code selber ansehen, ob noch was optimiert werden kann, kann auch 
helfen, oder von mir aus auch so ein Peep-Hole-Optimizer. Je nachdem, 
ein wenig ausprobieren kostet nicht viel.
Man kann ja auch selbst eine Statistik machen, welcher Code besser 
performt.

Lesen, also theoretische Überlegungen helfen natürlich auch noch, muss 
man halt auch lesen bzw. durcharbeiten.
https://www.ripublication.com/ijaer18/ijaerv13n12_107.pdf
https://www.researchgate.net/publication/338076963_Study_of_Median_Filters_and_Fuzzy_Rules_for_Filtering_of_Images_Corrupted_With_Mixed_Noise
https://www.researchgate.net/publication/261091013_Bidimensional_median_filter_for_parallel_computing_architectures

[Detlef _. (detlef_a)]

Hi,

meines Kenntnisstands nach heißt 9 sample median: 9 samples der Größe 
nach sortieren und den Fünftgrößten nehmen. Wenn ein neuer kommt, den 
ältesten aus der Liste rausnehmen und den neuen einsortieren. Sortieren 
ist O(n*log(n)), Einsortieren in eine vorsortierte Liste ist vielleicht 
sogar nur O(log(n)). Das geht schnell.

Interessantes Problem, wenns morgen regnet probier ich das vllt. mal 
aus.

Cheers
Detlef

[Peter M. (r2d3)]

Hallo Detlef_.,

ob das denn wohl so sinnvoll ist, mit der O()-Notation für n=9 zu 
argumentieren? :)

[Detlef _. (detlef_a)]

Peter M. schrieb:
> Hallo Detlef_.,
>
> ob das denn wohl so sinnvoll ist, mit der O()-Notation für n=9 zu
> argumentieren? :)

Ja ebent, das geht immer schnell.

[Detlef _. (detlef_a)]

Hier hab ich einen mittelmäßig effizienten medianfilter gebaut, Laufzeit 
ist O(n). Es wird eine sortierte Liste mitgeführt. Es wird aber noch so 
einiges hin und hergeschoben, das ist unschön und kostet Zeit. Mit mehr 
Überlegen und mehr Listen könnte man die Laufzeit vllt. auf O(log(n)) 
drücken, aber, wie Peter ja richtig sagte, n ist sowieso klein.

Cheers
Detlef

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

/***************************************/
float median(float neu){
/***************************************/
int32_t k,m;

#define NN (7)

static float alt[NN];
static  int32_t gross[NN];
static int32_t w=0;

//Ini

if(w==0){
 for(m=0;m<NN;m++)alt[m]=0.0;
 for(m=0;m<NN;m++)gross[m]=m;
 w=1;
}

// Den ältesten rausschmeissen
k=0;
while(1) {
 if(gross[k]==NN-1) break;
 k++;
}

//printf("alt %d  ",k);

while(1) {
 if(k==0) break;
 gross[k]=gross[k-1];
 k--;
}

// Wo musser denn hin?
k=0;
while(1) {
 if(alt[gross[k+1]]>neu) break;
 k++;
 if(k==NN-1) {break;};
}
//printf("hoch %d      ",k);

// Platz da
m=0;
while(1) {
 if(k==m) break;
 gross[m]=gross[m+1];
 m++;
}

gross[k]=-1;
//printf("gr %d %d %d                   ",gross[0],gross[1],gross[2]);
for(m=0;m<NN;m++) gross[m]++;
//printf("gr %d %d %d                   ",gross[0],gross[1],gross[2]);

// ruecken
m=NN-1;
while(1) {
 alt[m]=alt[m-1];
 m--;
 if(m==0) break;
}

alt[0]=neu;

return(alt[gross[(NN-1)/2]]);
}


/***************************************/
int main()
/***************************************/
{ int32_t k;
  float bla;
  float ee;

    for(k=0;k<20;k++){
      bla=k+1;
      if(k==13) bla=1000;
      ee=median(bla);
      printf("%d %.0f  \r\n",k,ee);
    }
    return 0;
}

[Christoph M. (mchris)]

>Hier hab ich einen mittelmäßig effizienten medianfilter gebaut, Laufzeit
>ist O(n).

Danke für den Code. Wenn ich darauf schaue, müsste die Laufzeit gegen

$$\infty$$

 gehen.


Das scheint mir eher eine Debug-Version zu sein, im Code steht 
haufenweise:

while(1)

[Christoph M. (mchris)]

Hier mal ein Benchmark mit Float.

// install median library:
// https://github.com/RobTillaart/RunningMedian

#include <RunningMedian.h>

RunningMedian samples = RunningMedian(3);

typedef float median_t;

void setup()
{
  Serial.begin(115200);
}

void loop()
{
  uint32_t ts;
  median_t x = random();
  median_t y;

  ts = micros();
  //y = median(x);
  samples.add(x);
  y = samples.getMedian();  
  uint32_t ende = micros();
  
  Serial.println("value" + String(y) + " time used [us]: " + String(ende - ts));
  delay(10);

}

[Christoph M. (mchris)]

Zum Vergleich die leere Medianfunktion.

[Norbert (der_norbert)]

Christoph M. schrieb:
> Das scheint mir eher eine Debug-Version zu sein, im Code steht
> haufenweise:
> while(1)

Vielleicht machst du dich mal mit ›break‹ vertraut.  ;-)

[Klaus (feelfree)]

Detlef _. schrieb:
> // Wo musser denn hin?
> k=0;
> while(1) {
>  if(alt[gross[k+1]]>neu) break;
>  k++;
>  if(k==NN-1) {break;};
> }

Süß. Meine erstes Bubblesort auf dem C64 vor genau 40 Jahren war ähnlich 
gut lesbar.

[Christoph M. (mchris)]

>Vielleicht machst du dich mal mit ›break‹ vertraut.  ;-)

Upps, ich bin damit vertraut, aber trotz angestrengtem Hinschauen hat es 
mein Hirn ignoriert.

Die Detlev'sche Version ist etwas schneller als die Tillaart'sche :-)
Ob der Median richtig sortiert, habe ich nicht geprüft, nur wie schnell 
der Algo läuft ..

// Algo from Detlef _. (detlef_a) 19.05.2024 20:53
// https://www.mikrocontroller.net/topic/567716#7668204

typedef float median_t;

#define NN (3)

median_t median(median_t neu)
{
  int32_t k, m;
  static median_t alt[NN];
  static int32_t gross[NN];
  static int32_t w = 0;

  //Ini
  if (w == 0) 
  {
    for (m = 0; m < NN; m++) alt[m] = 0.0;
    for (m = 0; m < NN; m++) gross[m] = m;
    w = 1;
  }

  // Den ältesten rausschmeissen
  k = 0;
  while (1) 
  {
    if (gross[k] == NN - 1) break;
    k++;
  }

  while (1) 
  {
    if (k == 0) break;
    gross[k] = gross[k - 1];
    k--;
  }

  // Wo musser denn hin?
  k = 0;
  while (1) 
  {
    if (alt[gross[k + 1]] > neu) break;
    k++;
    if (k == NN - 1)      break;
  }

  // Platz da
  m = 0;

  while (1) 
  {
    if (k == m) break;
    gross[m] = gross[m + 1];
    m++;
  }

  gross[k] = -1;
  for (m = 0; m < NN; m++) gross[m]++;
  // ruecken
  m = NN - 1;
  while (1) 
  {
    alt[m] = alt[m - 1];
    m--;
    if (m == 0) break;
  }

  alt[0] = neu;
  return (alt[gross[(NN - 1) / 2]]);
}

void setup()
{
  Serial.begin(115200);
}

void loop()
{
  uint32_t ts;
  median_t x = random();
  median_t y;
  
  ts = micros();
  y = median(x);
  uint32_t ende = micros();
  
  Serial.println("value" + String(y) + " time used [us]: " + String(ende - ts));
  delay(10);

}

[Norbert (der_norbert)]

Hab's gerade mal in – Obacht – Python programmiert und auf dem 
Mikrocontroller laufen lassen.

FILTERSIZE = 9
same values:       10.21 µs
linear up:         24.76 µs
linear down:       23.29 µs
random values:     24.09 µs

Gar nicht mal so übel. Kommt zwar nicht an die Thumb Assembler Version 
heran, aber immerhin gut für ~40.000 Werte pro Sekunde…

[Klaus (feelfree)]

Norbert schrieb:
> FILTERSIZE = 9

Christoph M. schrieb:
> #define NN (3)

Äpfel und Birnen?

[Norbert (der_norbert)]

Klaus schrieb:
> Norbert schrieb:
>> FILTERSIZE = 9
>
> Christoph M. schrieb:
>> #define NN (3)
>
> Äpfel und Birnen?

Nein Klaus, viel einfacher. Einmal Drei und einmal Neun. Kann man auch 
mit etwas Eigeninitiative ablesen.

Hint: Wenn man's mit ›3‹ macht, dann wird's schneller.
Wenn man's mit ›15‹ macht, dann wird's langsamer.

Hint2: Es ging nicht um einen direkten Vergleich.

Hint3: Bei einer Filterlänge von ›3‹ wäre der Aufwand des beschriebenen 
Algorithmus um Größenordnungen zu hoch!
Das könnte man mit einigen wenigen simplen Vergleichen viel schneller 
haben.

[Klaus (feelfree)]

Christoph M. schrieb:
> Upps, ich bin damit vertraut, aber trotz angestrengtem Hinschauen hat es
> mein Hirn ignoriert.

Ist auch ein bisschen ungewöhnlich bis pervers, die Abbruchbedingung in 
ein if zu packen, statt direkt in das while....

[Christoph M. (mchris)]

>Hab's gerade mal in – Obacht – Python programmiert und auf dem
> Mikrocontroller laufen lassen.

Könntest du den Python-Code posten? Ich würde es gerne "nacheruieren".
Welchen Mikrocontroller hast du verwendet? Da gibt es große Unterschiede 
in der Geschwindigkeit.

[Christoph M. (mchris)]

Hier mal die Laufzeit auf dem Atmega328 mit int16-Werten und Indexlänge 
uint8.
Zieht man die Laufzeit der leeren Funktion von oben (~4us) ab, kommt man 
auf eine durchschnittliche Laufzeit von 6us.
Nicht schlecht, würde ich sagen.

typedef int16_t median_t;
typedef uint8_t index_t;

#define NN (3)

median_t median(median_t neu)
{
  index_t k, m;
  static index_t gross[NN];
  static index_t w = 0;
  static median_t alt[NN];

[Klaus (feelfree)]

Christoph M. schrieb:
> #define NN (3)

Dafür braucht es keine median-Funktion mit sage und schreibe 5 
while-Schleifen.

[Christoph M. (mchris)]

>Dafür braucht es keine median-Funktion mit sage und schreibe 5
>while-Schleifen.

Jetzt besser?:

#define NN (11)

[Klaus (feelfree)]

Christoph M. schrieb:
> Jetzt besser?:

Genau.

Vielleicht kannst Die Laufzeit mal mit einer einfachen straight-forward 
Implementierung vergleichen?

Ich hätte das so implementiert:

int cmpfunc (const void* a, const void* b) {
   return ( *(median_t*)a - *(median_t*)b );
}

median_t median(median_t neu)
{
  static int32_t idx = 0;
  static median_t timebased[NN];
  static median_t sorted[NN];
  
  timebased[idx++] = neu;
  idx %= NN;

  memcpy(sorted, timebased, sizeof(sorted));
  qsort(sorted, NN, sizeof(median_t), cmpfunc);

  return (sorted[(NN - 1) / 2]);
}

[Christoph M. (mchris)]

Klaus (feelfree)
>Vielleicht kannst Die Laufzeit mal mit einer einfachen straight-forward
>Implementierung vergleichen?

Deine Implementierung ist auf jeden Fall schön übersichtlich. Ich habe 
ein wenig gezögert, weil ich vermutete, dass die in der in der Arduino 
IDE verwendete Compiler keine STD-Lib und damit kein QSort kennt. Aber 
siehe da, es kompiliert out of the Box.
Die Performance ist allerdings nicht so gut wie bei der Detlev'schen 
Version.
Hiermit möchte ich Detlev für seinen Beitrag ausdrücklich loben.

[Detlef _. (detlef_a)]

Hallo, vielen Dank. Ich bin noch nicht zufrieden. Dieses Umschaufeln ist 
nicht cool. In eine sortierte Liste ein Element neu einsortieren geht 
mit 'linear merge' in O(log log n), mein hochbubbeln ist O(n).Das 
älteste Element fliegt raus. Ich überlege an einer Datenstruktur die 
einerseits immer direkten Zugriff auf das älteste Element gibt aber 
andererseits schnelles Einsortieren des Neuen erlaubt. Wenn der Neue 
kommt werden alle eins älter, alle Altersangaben eins hochzählen ist 
aber O(n) .
Hm
Cheers
Detlef

[Klaus (feelfree)]

Detlef _. schrieb:
> Wenn der Neue
> kommt werden alle eins älter, alle Altersangaben eins hochzählen ist
> aber O(n)

Das Alter durch stures inkrementieren darstellen? Ein neues Element 
bekommt eine höhere Zahl. der älteste hat dann die kleinste Zahl.

[Christoph M. (mchris)]

Vielleicht könnte eine "Linked List" helfen, auf die zusätzlich ein 
Pointer mit dem Alter zeigt. Ich vermute aber, dass bei kurzen 
Medianfiltern dann eher der Verwaltungsaufwand zu groß ist.

[Klaus (feelfree)]

Hätte nicht gedacht, dass qsort() bei kurzen Listen soo schlecht 
performt....

Hier mal eine Version ohne qsort.
Eine binäre Suche nach der einzusortierenden Wert lohnt sich glaube ich 
bei so kurzen Listen kaum.
Bleibt immer noch die unschöne lineare Suche nach dem ältesten Wert.

median_t median3(median_t new_value)
{
  static uint8_t initialized = 0;
  static int32_t act_age = 0;

  struct listelement
  {
    median_t value;
    uint8_t age;
  };
  static struct listelement list[NN];

  uint8_t i;
  uint8_t index_new;
  uint8_t index_to_delete;

  if (!initialized)
  {
    for (i=0; i<NN; i++)
       list[i].age = i;
    initialized = 1;
  }

  for (i=0; i<NN; i++)
  {
    if (new_value < list[i].value) break;
  }
  index_new = i;

  for (i=0; i<NN; i++)
  {
    if (list[i].age == act_age) break;
  }
  index_to_delete = i;

  if (index_to_delete < index_new)
  {
    index_new--;
    // move elements [index_to_delete+1 .. index_new] to the left
    for (int8_t idx=index_to_delete+1; idx<=index_new; idx++)
        list[idx-1] = list[idx];
  }
  else if (index_new < index_to_delete)
  {
    // move elements [index_new .. index_to_delete-1] to the right
    for (int8_t idx=index_to_delete-1; idx>=index_new; idx--) 
        list[idx+1] = list[idx];
  }

  list[index_new].value = new_value;
  list[index_new].age = act_age;

  act_age++;
  act_age %= NN;

  return (list[(NN - 1) / 2].value);
}