STM32F405 - DSP Medianfilter

[Lea R. (michael1909)]

Hallo,

hat von euch schon mal jemand einen Medianfilter auf einen STM32F4 
programmiert?

Ich benötige diesen um aus einer Datenreihe Ausreißer zu entfernen.

Meinen derzeitigen Medianfilter habe ich mit for-schleifen gelöst.
mit einer Filterlänge von 5 Punkte.

kennt ihr eine DSP-Funktion des Cortex-M4 einen Median-Filter zu 
optimieren?

Mein Ziel ist es, eine bessere Performance zu erreichen.
Derzeit benötigt mein Medianfilter bei einer Datenlänge von ca. 
1000Punkte ungefähr 6ms.

LG Michael

[Gerald G. (gerald_g)]

Michael Höller schrieb:
> Meinen derzeitigen Medianfilter habe ich mit for-schleifen gelöst.

Dass da irgendwo ne Schleife ist, sollte klar sein.
Wenn etwas optimiert werden soll, brauchen wir den code. Ansonsten 
findet man sogar auf dieser Seite Implementierungen in c

[Detlef K. (adenin)]

Wenn man den Ist-Zustand nicht kennt, dann kann man keine brauchbaren 
Tipps zur Optimierung geben.

Beispiel für völlig unbrachbaren Tipp: :)
Du schreibst
>Derzeit benötigt mein Medianfilter bei einer Datenlänge von ca.
>1000Punkte ungefähr 6ms.

Dann schreib ich einfach: Verzehnfache deine Taktfrequenz, dann brauchst 
Du nur noch 0.6ms.

[Lea R. (michael1909)]

:) Okay!
hier mein derzeitiger Programmcode:


static void MedianFilter(void)
{
    // TODO
    //

    static unsigned short z = 0, y = 0, x = 0, v = 0, min = 0, temp = 0;
    static unsigned short window[5];


    for(z = 2; z < (buffer_index - 2); z++)
    {
        for(y = 0; y < 5; y++)
        {
            window[y] = FilterData.D_Buf_CACHE[y+z];
        }

        for(x = 0; x < 3; x++)
        {
            min = x;
            for((v = x +1); v < 5; v++)
            {
                if(window[v] < window[min])
                    min = v;
            }
            temp = window[x];
            window[x] = window[min];
            window[min] = temp;
        }
        // get result --> middle element
        FilterData.D_MF_Buf[z - 2] = window[2];
    }
}

Fertige DSP - Funktionen gibt es ja anscheinend nicht, ähnliche wie bei 
den FIR-Filtern zum beispiel!

vielen Dank

[Erich (Gast)]

Auf der Seite
http://stackoverflow.com/questions/480960/code-to-calculate-median-of-five-in-c-sharp
findet sich einiges genau dazu.
Die Antwort von
    answered Jan 26 '09 at 21:20
erscheint nachvollziehbar, während das von
    answered Jan 22 '10 at 12:00
doch ziemlich heftig ist.
Gruss

[Arc N. (arc)]

Erster Gedanke: Muss der Filter nachträglich laufen oder ginge das auch 
direkt beim Eintreffen von genügend (Fenstergröße) Daten?
Zweiter Gedanke: Das Sortieren ist in dieser Form nicht nötig. Stichwort 
ist find-k-smallest-element / Blum-Floyd-Pratt-Rivest-Tarjan-Algorithmus 
1) was den Aufwand an dieser Stelle auf O(n) senken würde (bei dieser 
Fenstergröße eigentlich nicht nötig).
Dritter Gedanke: Warum braucht der jetzige Algorithmus so lange?
Grob überschlagen bei 1000 Elementen:
- 1000 * 5 Lesezugriffe auf das Ausgangsarray, 5000
- 1000 * 5 Schreibzugriffe auf das Fenster, 5000
- 1000 * (18 Schreib- und 6 Lesezugriffe und 9 Vergleiche), 33000
- 1000 * (1 Schreib- und Lesezugriff zum Zurückschreiben), 2000
Grob 45000 Operationen in 6 ms also 7.5 MHz...

1) http://people.csail.mit.edu/rivest/pubs/BFPRT73.pdf
(wer mal Info oder ähnliches studiert (hat), wird wahrscheinlich mehr 
als einen der Autoren wiedererkennen ;))

[Gerald G. (Gast)]

Ich werfe mal etwas Wahrscheinlichkeit mit in die Auswertung.
Die Auswertung dauert also ca. 33 Cycles.
Bei

for(z = 2; z < (buffer_index - 2); z++)
    {
        for(y = 0; y < 5; y++)
        {
            window[y] = FilterData.D_Buf_CACHE[y+z];
        }

wird ja im Endeffekt Daten[y] bis Daten[y+4] untersucht. Beim nächsten 
Durchlauf untersucht man Daten[y+1] bis Daten[y+5]. Es wird also nur 
Datenpunkt Data[y] mit Data[y+5] ersetzt, der Rest bleibt erhalten.
Ist also der Datenpunkt Data[y] der herausfällt, nicht der Median, gibt 
es zwei Möglichkeiten: Er ist entweder größer als der Median oder 
kleiner. Trifft auf Data[y+5] das selbe zu (Also beide sind 
größer/kleiner als der Median) ändert sich der Median nicht.
Also:

 if((Data[y] > Median_alt) && (Data[y+5] > Median_alt))
       Median_neu = Median_alt;

Das selbe gilt auch für Data[y] < Median_alt

Fügt man also zusätzliche 4 Vergleiche zu den 33 Cycles der Auswertung, 
hat man 37 Cycles. Zu 40% (Bei statistischen Zahlen) trifft aber der 
Fall ein, dass Data[y] und Data[y+5] im selben Verhältnis zum Median 
stehen, wodurch sich der Durchlauf für diese Daten auf 4 Cycles 
reduziert.

Das ergibt im Durchschnitt:
0,4* 4 Cycles + 0,6 * 37 = 23,8 Cycles.

Man reduziert also die 45000 Cycles auf 35800 Cycles, also von 6ms auf 
4,77ms

Das wäre so meine Idee, ohne jetzt groß den Auswertealgorhithmus zu 
ändern :D

[Arc N. (arc)]

"Fieser" Trick:
Wenn die ersten 5 Werte sortiert sind, braucht bei den nächsten 5 Werten 
nur noch geschaut werden wo sich der neue einsortiert.

  for (int i = 1; i < 1000 - 5; i++) {
    // value to be deleted from window
    unsigned short vd = array[i - 1];
    // value to be added to window
    unsigned short va = array[i + 4];
    for (int j = 0; j < 5; j++) {
      if (window[j] == vd) {
        // total operations in k-loop: 5 + 5 comparisons + 4 read and 4 + 1 write operation
        for (int k = 0; k < 5; k++) {
          if (k == j) continue;
          if (va < window[k]) {
            // now insert va...
            //  d, v1, v2, v3, v4
            // v0,  d, v2, v3, v4
            // v0, v1,  d, v3, v4
            // v0, v1, v2,  d, v4
            // v0, v1, v2, v3, d            
          }
        }
        break;
      }
    }
  }

Das sollten, falls der Algorithmus so stimmt, dann 1000 * (5 + 5 + 5 + 4 
+ 5) = 1000 * 24 = 24000 Operationen sein
(da fehlen noch die Vergleiche, um die passenden Verschiebungen im 
Fenster zu erledigen)

[Lea R. (michael1909)]

vielen Dank für Eure Denkanstöße,

werde jetzt einmal ein paar Vorschläge testen!

Danke!!!
Michael.