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
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.
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
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.
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
Sowas macht man nicht in Fliesskomma. Dafür gibt es Festkommaformate wie Q15.1 etc.
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 ;)
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
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
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?
@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.
Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.