Digitaler Filter - Bandpass

[Dominik (Gast)]

Hallo,

ich möchte Neuland betreten und etwas über digitale Filter lernen. 
Konkret geht es mir darum, die Amplitude eines 50Hz-Signals zu 
bestimmen. Problem dabei: Das Signal ist nicht "sauber". Es enthält 
Oberwellen und schlimmstenfalls eine Überlagerung mit einem starken 
16,7Hz Signal von einem Oberleitungsfeld der Bahn.
Deshalb bräuchte ich normalerweise wohl einen (analogen) Bandpass, um 
ein sauberes 50Hz Signal wieder daraus zu filtern. Sowas soll aber auch 
digital gehen. Nur habe ich ehrlich gesagt keine Vorstellung, wie sowas 
funktioniert. Ziemlich ratlos lese ich den Wikipediaartikel: 
http://de.wikipedia.org/wiki/Digitales_Filter
Hat vielleicht jemand einen Link mit einer einfachen Einführung in die 
Thematik? Ich möchte die 50Hz-Amplitude über ADC-Wandlungen eines 
einfachen Mikrocontroller (AtMega oder PIC16) auswerten. Wäre so ein 
Mikrocontroller programmiertechnisch auch in der Lage, dabei gleich 
einen digitalen 50Hz Bandpass anzuwenden? Gibt es vielleicht 
Beispielprojekte, wo ein digitaler Bandpass (Frequenz sei mal egal) mit 
einem Mikrocontroller realisiert ist?

Danke im voraus,
Dominik

[Lars S (Gast)]

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.

[Lars S (Gast)]

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

[Beobachter (Gast)]

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

[OR (Gast)]

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

[Helmut L. (helmi1)]

Dominik schrieb:
> Ich möchte die 50Hz-Amplitude über ADC-Wandlungen eines
> einfachen Mikrocontroller (AtMega oder PIC16) auswerten. Wäre so ein
> Mikrocontroller programmiertechnisch auch in der Lage, dabei gleich
> einen digitalen 50Hz Bandpass anzuwenden?

Klar geht das. Mit einem ATMega88 mit 16MHz Takt kann man einen 
digitalen IIR Bandpass bis 2KHz bauen (5 KHz Samplerate). Einen 50Hz 
Bandpass mit Guetefaktor von 1000 ist da kein Problem.

[Armin (Gast)]

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!

[Dominik (Gast)]

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

[Helmut L. (helmi1)]

Sowas macht man nicht in Fliesskomma. Dafür gibt es Festkommaformate wie 
Q15.1 etc.

[Lars S (Gast)]

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;

[Lars S (Gast)]

>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

[Lars S (Gast)]

so etwa (nicht ausprobiert)

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

int sinusanteil;
int cosinusanteil;
int sinfiltbuffer[20];
int cosfiltbuffer[20];
long sinfilt = 0;
long cosfilt = 0;
int ergebnis_im_quadrat;

Task1ms () {
  static int i = 0;

  // Phasenzähler von 0..19
  i++;
  if (i==20)
    i=20;

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

  // Gleitender Mittelwert Sinusanteil
  sinfilt -= sinfiltbuffer[i];
  sinfiltbuffer[i] = sinusanteil;
  sinfilt += sinfiltbuffer[i];

  // Gleitender Mittelwert Cosinusanteil
  cosfilt -= cosfiltbuffer[i];
  cosfiltbuffer[i] = cosinusanteil;
  cosfilt += cosfiltbuffer[i];

  // Geometrischer Mittelwert aus sin und cos
  ergebnis_im_quadrat = sinfilt/20*sinfilt/20 + cosfilt/20*cosfilt/20;

}

[Dominik (Gast)]

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 ;)

[Marius W. (mw1987)]

Schau dir mal den Goertzel-Algorithmus an. Das ist ein Algorithmus zur 
Fourier-Transformation einer einzigen Frequenz. Da das ganze lässt sich 
in Form eines rekursiven Filters implementieren.

Weitere Infos hier: http://de.wikipedia.org/wiki/Goertzel-Algorithmus

MfG
Marius

[Achim M. (minifloat)]

Kleiner Ausflug in die Systemtheorie:

Ich habe ein Signal x(t), das schicke ich durch ein System mit der 
Impulsantwort g(t) und erhalte das Ausgangssignal y(t).
Die mathematische Formulierung ist:
y(t) = x(t) * g(t)
Mit * meine ich den Faltungsoperator.
In Wirklichkeit wird also gefaltet, das ist ein ziemlich krumm 
anmutendes Integral.

Besser gehts schon im Laplace-Bildbereich. Dort ist eine Faltung im 
Zeitbereich eine Multiplikation. Man muss mit einer Korrespondenztabelle 
seine Signale transformieren, kann sie dann multiplizieren und wieder 
zurücktransformieren. Y(s) = G(s) mal X(s).
Wenn man seine Übertragungsfunktion G(s) in einen Zähler Z(s) und einen 
Nenner N(s) aufspalten kann, kann man direkt ein Blockschaltbild 
angeben, welches nur Integratoren zur Berechnung des Ausgangssignales 
aus dem Eingangssignal benutzt. Dieses Blockschaltbild
gilt selbstverständlich auch für den Zeitbereich, da die 
Laplace-Transformation linear ist. Der Integrator ist der 
Elementarbaustein für zeitkontinuierliche Signalprozessoren, welcher 
sich aber bereits einfachst in Software bauen lässt: ein additiv zu 
verändernder, direkt auslesbarer Zwischenspeicher. Man muss nur die 
Integrierzeit mit der Abtastfrequenz
in Relation setzen und vor dem Addieren damit multiplizieren.

Oder eben vom Laplace-Bildbereich in den z-Bildbereich (zeitdiskreter 
Bildbereich) transformieren.
Dazu gibt es mehrere Transformationen, die je nach Signalcharakter 
unterschiedlich genau im Ergebnis sind. Es gibt die Sprunginvariante, 
die Impulsinvariante, die Bilineare Transformation, Vorwärts- und 
Rückwärtsdifferentiation...
Auch hier kann man das Blockschaltbild nach Auflösung in Z(z) und N(z) 
direkt in ein Blockschaltbild umwandeln. Dieses besitzt nur 
Verzögerungsglieder in Form von Totzeitglidern(Zwischenspeicher für 
Werte).
Die Totzeit ist gleich deiner Abtastzeit, also der Kehrwert deiner 
Abtastfrequenz. So ein Blockschaltbild lässt sich sehr einfach in 
Software gießen.

mfg mf

[Lars S (Gast)]

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:

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

int sinusanteil;
int cosinusanteil;
long sinfilt = 0;
long cosfilt = 0;
int ergebnis_im_quadrat;
int i=0, ii =0;

Task1ms () {

  // Phasenzähler von 0..19
  i++;
  if (i==20)
    i=0;

  // Filterzähler von 0..99ms
  ii++;

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

  // alle 100ms ein Ergebnis
  if (ii==100) {

    // Geometrischer Mittelwert aus sin und cos
    ergebnis_im_quadrat = (sinfilt*sinfilt + cosfilt*cosfilt) / 100;
    
    sinfilt =0;  //Mittelwertbildung wieder zurücksetzen
    cosfilt =0;
  }
  else {
   // Mittelwert bilden
    sinfilt += sinusanteil;
    cosfilt += cosanteil;
  }


}

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

sinfilt_long += sinusanteil - sinfilt;
sinfilt = sinfilt_long / 1024;

[Lars S (Gast)]

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)

[OR (Gast)]

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

[Lars S (Gast)]

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.

[Dominik (Gast)]

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?

[Helmut L. (helmi1)]

@Dominik

Falls du es mal mit einem Digitalen IIR Filter probieren moechtest ..

In Tonedecoder ist das IIR Filter enthalten. Zur berechnung der 
Fractionalmultiplikation habe ich die Multiplikationroutine in Assembler 
geschrieben Q22.s  (24 Bit x 24 Bit Multiplikation die es in C nicht 
gibt)
In C wuerde die zu lange dauern.

Das IIR Filter ist im Code auf 50Hz Mittenfrequenz und 1 Hz Bandbreite 
eingestellt.