Digitaler Filter - Bandpass

Erzeuge Dir (mathematisch) ein Sinus und Cosinussignal mit 50Hz.
Multipliziere das jeweils mit deinem ADC Wert.
Integriere das dann für n Sekunden auf. Bilde den Quadratischen 
Mittelwert aus beiden Ergebnissen.

Dann hast Du das Ergebnis, als ob Du ein Bandfilter mit Bandbreite 1/n 
hättest.

Achja. Das Ergebnis ist also sowas ähnliches wie bei ner FFT. Nur eben 
für exakt eine Frequenz.

ist meines Wissens aber sehr steil flankig.
Nennt sich DFT.

Die Steilflankigkeit hängt an der Fensterbewertung. Der vorgeschlagene 
Algorithmus ist ein Rechteckfenster. Mann kann aber nach der 
sincos-Multiplikation auch Tiefpässe mit geeignetem Frequenzgang wählen. 
Dadurch kann man fast beliebige Bandpässe bauen. Ergebnis ist immer ein 
real-und Imaginärteil. Das ganze kann man auch sehr schön im 
Matlab/Octave ausprobieren.

Anm.: Falls man am Betrag des Bandpassausganges interessiert ist, dann 
ist noch eine geometrische Addition der beiden Stränge notwendig.

Grüße
OR

du kannst jede beliebige kontinuierliche lineare Übertragungsfunktion 
nehmen, deren Filtereigenschaften über ein Bode-Diagramm abzulesen sind.
--> http://de.wikipedia.org/wiki/Bode-Diagramm

dürfte bekannt sein.


Diese Übertragungsfunktion lässt sich mittels Tusin-Transformation in 
der Zeit diskretisieren:
--> 
http://de.wikipedia.org/wiki/Bilineare_Transformation_(Signalverarbeitung)


wenn man das auflöst und umformt, wird es zu einem digitalen 
Alogrithmus, der sich wunderbar implementieren lässt, WENN der in 
gleichmäßigen Abständen ausgeführt wird. Also Timer verwenden!

Lars S schrieb:
> Erzeuge Dir (mathematisch) ein Sinus und Cosinussignal mit 50Hz.
> Multipliziere das jeweils mit deinem ADC Wert.
> Integriere das dann für n Sekunden auf. Bilde den Quadratischen
> Mittelwert aus beiden Ergebnissen.
>
> Dann hast Du das Ergebnis, als ob Du ein Bandfilter mit Bandbreite 1/n
> hättest.

Das ist ja schon ne Menge Rechnerei. Ob das einem einfachen 
Mikrocontroller packen kann? Erstmal brauche ich also eine Sinus und 
Cosinusfunktion, die allein wahrscheinlich schon einige Rechenzeit 
beantspruchen.
Angenommen ich mache mit einem Timer 1000 ADC-Wandlungen pro Sekunde 
(sprich Abtastfrequenz 1kHz), sprich 20 Wandlungen pro Periode und 
berechne entsprechend dann auch 1000 50Hz-Sinus und -Cosinus-Werte pro 
Sekunde. Kann ich mit der Berechnung zu einem willkürlichen Zeitpunkt 
anfangen oder muss das noch irgendwie mit dem Nutzsignal synchronisiert 
werden?
Ich möchte beispielsweise 1 Sekunde integerieren um eine Bandbreite 1/n 
= 1/1s = 1Hz (sprich Bandpass von 49HZ-51Hz) zu bekommen.

obiges Beispiel mal Schritt für Schritt betrachtet:

1. ADC-Wert: Multiplikation mit aktuellen Sinus SIN(t) und Cosinus 
COS(t) und speichere das Ergebnis als Integral
Integral_1SIN = 1ms * SIN(t) * ADC-Wert
Integral_1COS = 1ms * COS(t) * ADC-Wert

2. ADC-Wert: Wieder Multiplikation mit aktuellen Sinus und Cosinus und 
speichere
Integral_1SIN = Integral1SIN + 1ms * SIN(t) * ADC-Wert
Integral_1COS = Integral1COS + 1ms * COS(t) * ADC-Wert
Integral_2SIN = 1ms * SIN(t) * ADC-Wert
Integral_2COS = 1ms * COS(t) * ADC-Wert

3. ADC-Wert: Wieder Multiplikation mit aktuellen Sinus und Cosinus und 
speichere
Integral_1SIN = Integral1SIN + 1ms * SIN(t) * ADC-Wert
Integral_1COS = Integral1COS + 1ms * COS(t) * ADC-Wert
Integral_2SIN = Integral2SIN + 1ms * SIN(t) * ADC-Wert
Integral_2COS = Integral2COS + 1ms * COS(t) * ADC-Wert
Integral_3SIN = 1ms * SIN(t) * ADC-Wert
Integral_3COS = 1ms * COS(t) * ADC-Wert

So bekomme ich 2*1000 Fliesskomma-Variablen (logischerweise als Array) 
als Puffer. Ganz schön viel! Das wird der Speicher des Mikrocontrollers 
schon knapp. Vielleicht muss ich die Bandbreite doch etwas grösser 
wählen: 500Werte, sprich +-2Hz Bandbreite tun es vielleicht auch. Aber 
bleiben wir mal bei dem Beispiel:

1000. ADC-Wert: Jetzt sei eine Sekunde vergangen, sprich die Variablen 
Integral_1SIN und Integral_1COS haben die Summe von 1000 Werten. Nun 
bilde ich das quadratische Mittel aus beiden Integralen sprich:
WURZEL ( (Integral_1SIN * Integral_1SIN + Integral_1COS * Integral_1COS) 
/ 2 )

Und dieses Quadratische Mittel kann ich nun als Momentanwert des 
Eingangssignals, welches durch einen digitalen Bandpass mit der schmalen 
Bandbreite von 1Hz gefiltert wurde, betrachten?

Nachteil zum analogen Filter wäre dann wohl, dass ich das gefilterte 
Eingangssignal (=Nutzsignal) in dem Beispiel erst mit einer Sekunde 
Verzögerung zur Verfügung hätte?
Wäre meine Rechung soweit richtig oder habe ich etwas falsch verstanden?


Vielen Dank,
Dominik

Nein. Das geht sehr einfach. Du brauchst da fast keine Resourcen. Auch 
kann man alles mit integer/long rechnen.

Hab das als Beispiel mal in einer Excel-Tabelle gemacht.

Sinus und Cosinus legst Du in einer Tabelle ab. Da du mit 1kHz abtastest 
und 50Hz Sin und Cos erzeugen willst, brauchst Du also nur eine Tabelle 
mit 20 Einträgen.

Du rechnest:

Sinusanteil = SinusAusTabelle()*ADCWert;
Cosinusanteil = CosinusAusTabelle()*ADCWert;

Sinusgefiltert = GleitenderMittelwertderletzten20Sinusanteile();
Cosinusgefiltert = GleitenderMittelwertderletzten20Cossinusanteile();

QuadratdesErgebnisses = Sinusgefiltert*Sinusgefiltert + 
Cosinusgefiltert*Cosinusgefiltert;

>Nachteil zum analogen Filter wäre dann wohl, dass ich das gefilterte
>Eingangssignal (=Nutzsignal) in dem Beispiel erst mit einer Sekunde
>Verzögerung zur Verfügung hätte?

Hast Du beim analogen Filter auch.
Du hast immer eine Zeitverzögerung mit 1/Bandbreite

so etwa (nicht ausprobiert)

1

int sintab[20]={0,10,19,26,30,32,30,26,19,10,0,-10,-19,-26,-30,-32,-30,-26,-19,-10};

2

int costab[20]={32,30,26,19,10,0,-10,-19,-26,-30,-32,-30,-26,-19,-10,0,10,19,26,30};

3

4

int sinusanteil;

5

int cosinusanteil;

6

int sinfiltbuffer[20];

7

int cosfiltbuffer[20];

8

long sinfilt = 0;

9

long cosfilt = 0;

10

int ergebnis_im_quadrat;

11

12

Task1ms () {

13

  static int i = 0;

14

15

  // Phasenzähler von 0..19

16

  i++;

17

  if (i==20)

18

    i=20;

19

20

  // Sinus und Cosinusanteil

21

  sinusanteil = adcvalue * sintab[i] / 32;

22

  cosinusanteil = adcvalue * costab[i] / 32;

23

24

  // Gleitender Mittelwert Sinusanteil

25

  sinfilt -= sinfiltbuffer[i];

26

  sinfiltbuffer[i] = sinusanteil;

27

  sinfilt += sinfiltbuffer[i];

28

29

  // Gleitender Mittelwert Cosinusanteil

30

  cosfilt -= cosfiltbuffer[i];

31

  cosfiltbuffer[i] = cosinusanteil;

32

  cosfilt += cosfiltbuffer[i];

33

34

  // Geometrischer Mittelwert aus sin und cos

35

  ergebnis_im_quadrat = sinfilt/20*sinfilt/20 + cosfilt/20*cosfilt/20;

36

37

}

Erst mal vielen Dank für Deine Mühe Lars!

Nur leider kann ich es immer noch nicht nachvollziehen, was Du da 
eigentlich machst. Ich habe mich die letzte Stunde mal näher mit Deiner 
Excel-Tabelle beschäftigt: Als Ergebnis bekommt man ja immer nahezu 1 
raus. Nur was soll das Ergebnis (jene 1) überhaupt darstellen? Es ist ja 
gerade nicht der durch den Bandpass gefilterte Sinus, den ich eigentlich 
erwartet hätte. Soll das die Amplitude sein? Lege ich eine andere 
Frequenz statt des 50Hz Sinus als Eingang, so kriege ich ein sehr 
schwankendes Ergebnis als quadratisches (oder geometrisches? Das steht 
wohl fälschlicherweise in der Exceltabelle als Spaltenüberschrift?) 
Mittel, aber nicht annährend 0, wie ich es von so einem steilen Bandpass 
eigentlich erwartet hätte.

Diese Fourier Transformationen sind mir doch ein wenig zu hoch. Gibt es 
dazu vielleicht ein gutes Buch zur Einführung, dass auch für mich 
(Nicht-Akademiker - nur Hobby Elektroniker) verständlich sein könnte? 
Grundlagen der digitalen Signalverarbeitung sind mir schon bekannt 
(Fourier-Analyse, Shannon-Theorem,PCM). Aber Filter sind für mich 
absolutes Neuland.
Nichtsdestotrotz: Einen fertigen Algorithmus, wie Du ihn beschreibst, 
kriege ich programmtechnisch gesehen vielleicht noch irgendwie umgesetzt 
- selbst wenn ich noch nicht verstehe, was der eigentlich machen soll ;)

Hallo Dominik,
das mit dem gleitenden Mittelwert war nicht ganz korrekt.
Anbei nochmals die Berechnung in Calc/Excel.

Ich hab das so erweitert, dass man jetzt auch die 16.6Hz als Störung 
hinzufügen kann.

Filterung ist nicht mehr ein gleitender Mittelwert, sondern entweder der 
Mittelwert nach 100ms (also Mittelwert aus 100 Messungen) (1. grüne 
Zahl), nach einer Sekunde (2. grüne Zahl), oder der gefilterte Wert mit 
einem PT1 Filter (Diagramm rechts unten).

Die Umsetzung als Programm wird dadurch auch nochmals einfacher. 
Beispiel:

1

int sintab[20]={0,10,19,26,30,32,30,26,19,10,0,-10,-19,-26,-30,-32,-30,-26,-19,-10};

2

int costab[20]={32,30,26,19,10,0,-10,-19,-26,-30,-32,-30,-26,-19,-10,0,10,19,26,30};

3

4

int sinusanteil;

5

int cosinusanteil;

6

long sinfilt = 0;

7

long cosfilt = 0;

8

int ergebnis_im_quadrat;

9

int i=0, ii =0;

10

11

Task1ms () {

12

13

  // Phasenzähler von 0..19

14

  i++;

15

  if (i==20)

16

    i=0;

17

18

  // Filterzähler von 0..99ms

19

  ii++;

20

21

  // Sinus und Cosinusanteil

22

  sinusanteil = adcvalue * sintab[i] / 32;

23

  cosinusanteil = adcvalue * costab[i] / 32;

24

25

  // alle 100ms ein Ergebnis

26

  if (ii==100) {

27

28

    // Geometrischer Mittelwert aus sin und cos

29

    ergebnis_im_quadrat = (sinfilt*sinfilt + cosfilt*cosfilt) / 100;

30

31

    sinfilt =0;  //Mittelwertbildung wieder zurücksetzen

32

    cosfilt =0;

33

  }

34

  else {

35

   // Mittelwert bilden

36

    sinfilt += sinusanteil;

37

    cosfilt += cosanteil;

38

  }

39

40

41

}

oder als PT1 Filter mit Tau = 1.024 Sec, also Bandpass mit 1Hz:

1

sinfilt_long += sinusanteil - sinfilt;

2

sinfilt = sinfilt_long / 1024;

Ups. Tabelle vergessen.

Dreh mal ganz rechts oben an den beiden Werten "Rauschanteil" und 
"Anteil 16.6Hz" und schau, wie sich das Ergebnis ( die beiden grünen 
Zahlen verändern)

Hallo,

die gleitende Mittelwertbildung halte ich für etwas zu aufwendig.
prinzipiell stelle ich mir das so vor:

// Sinus und Cosinusanteil
  sinusanteil = adcvalue * sintab[i] ;
  cosinusanteil = adcvalue * costab[i] ;

// glaetten mit Hilfe eines PT1
   sinusanteil_g += (sinusanteil-sinusanteil_g)>>4; // durch 16 teilen
   cosinusanteil_g += (cosinusanteil-cosinusanteil_g)>>4;

// quadratischen Bandpassausgang, kann man auch in einer langsameren 
Zeitscheibe machen:
   res = sinusanteil_g * sinusanteil_g + cosinusanteil_g*cosinusanteil_g

Grüße
OR

Hallo OR,

ja - es ist richtig, dass der gleitende Mittelwert zu aufwendig ist. Das 
sehe ich mittlerweile genauso.

Der Algorithmus des PT1 den Du beschreibst hat allerdings immer eine 
Abweichung vom Endwert von 16 zur Folge. Der 2 Postings darüber 
beschrieben PT1 Algorithmus hat das nicht. Dafür benötigt man eine 
zusätzliche long-Variable.

Muss man abwägen, was man will.

Erst mal vielen vielen Dank! Die Algorithmen funktionieren jetzt bestens 
und sind echt einfacher, als ich gedacht hätte!
Das beschriebene Verfahren ist aber nur bei genau einer Frequenz 
irgendwie sinnvoll, indem es dessen Amplitude bestimmt (was genau auf 
mein Ausgangsproblem zutrifft!). Insofern hat es aber mit einem 
"digitalen Bandpass" noch wenig zu tun.
Wenn ich (nur rein theoretisch) einen breiteren Frequenzbereich als 
Bandpass haben möchte (z.B. 100-300Hz), worin mich besonders 
Frequenzschwankungen interessieren, dann greift das Verfahren leider 
nicht mehr. Da interessieren mich ja nicht die Amplituden (die bei 
Tonübertragungen oder ähnlichem zudem eh sehr schwankend sind), sondern 
die Momentanwerte. Ginge so etwas auch in ähnlicher Weise, so dass ich 
als Ergebnis Momentanwerte erhalte?