Forum: Digitale Signalverarbeitung / DSP / Machine Learning Besselfilter umsetzen

http://www.mikrocontroller.net/articles/Digitalfilter_mit_ATmega

Hallo,

vielen dank für die Antwort, aber irgendwie habe ich es immer noch net 
verstanden, wie ich das umsetzten soll.

Es muss doch eine allg. gültige Gleichung geben, bei der die Abtast oder 
Grenzfrequez mit einfließen. Ist doch in der Analog Technik nicht 
anders.

Kann mir jemand einen Tipp geben?

Gruß

Hallo Ruslan,

die Ermittlung der Koeffizienten eines Besselfilters ist nicht trivial, 
daher verwendet man entweder dedizierte Designroutinen (Matlab?) oder 
schlägt sie in Tabellen nach (z.B. Halbleiterschaltungstechnik von 
Tietze/Schenk). Dort findet man Dimensionierungshilfe für die analoge 
Implementierung, die in Differenzengleichungen überführt werden will.

Warum ist das so kompliziert? Mit Besselfiltern versucht man, den 
linearen Phasengang, also die konstante Gruppenlaufzeit von FIR-Filtern 
mit IIR-Filtern zu realisieren. Es handelt sich jedoch nur um eine 
Approximation der konstanten Gruppenlaufzeit. Zur Ermittlung der 
Koeffizienten für höhere Ordnungen ist ein nichtlineares 
Gleichungssystem zu lösen. Schließlich gibt es Rekursionsformeln zur 
Koeffizientenberechnung, die die Bessel-Polynome enthalten; daher der 
Name.

Woher kommt Deine Anforderung, unbedingt ein Bessel-Filter 2. Ordnung 
mit 400Hz Grenzfrequenz einzusetzen?
Benötigst Du wirklich eine konstante Gruppenlaufzeit?
Hast Du harte Echtzeitanforderungen, bei denen einige zig Abtastwerte 
Verzögerung stören (kritische Totzeit in Regelung)?

Wenn ich z.B. Sensorwerte filtern möchte, setze ich meist ein ganz 
einfaches IIR-Filter erster Ordnung wie folgt ein:

y(n)= a * x(n) + (1-a)*y(n-1)

x(n) sind die Eingangswerte (bei Dir Beschleunigungswerte),
y(n) sind die gefilterten Ausgangswerte,
a = 1-exp(-2 pi fg/fs) mit fg als 3dB-Grenzfrequenz und fs als 
Abtastfrequenz.

Falls Gruppenlaufzeitverzerrungen wirklich stören, kannst Du auch ein 
symmetrisches FIR-Filter einsetzen. Die Anzahl der Koeffizienten und 
damit die konstante Gruppenlaufzeit (=Koeffizienten/fs/2) hängt von der 
nötigen Flankensteilheit ab. Die Gruppenlaufzeit im Durchlaßbereich des 
Filters wird dann höher als im IIR-Fall sein und könnte in einer 
Regelung durch die zu große Totzeit zu Problemen führen. Dafür ist das 
Design von FIR-Filtern sehr einfach.

Ich hoffe das hilft etwas.

Viele Grüße!

Gerrit, DL9GFA

Dazu gibt es Matlab und ähnliche Programme. Im kostenlosen Scilab zum 
Beispiel mit "eqiir"
http://www.scilab.org/product/man/eqiir.html
da gibts leider kein Bessel, nur Butterworth, Chebysheff und Cauer 
(elliptic), da wäre "butt" das nächstliegende.

Analog zum Notchfilter im Artikel müßte der Tiefpass in Scilab so 
lauten:

1

// Lowpass-Filter für 400 Hz, Samplerate =10000Hz

2

// Eckfrequenzen 400Hz, 800Hz:

3

omega=[2*%pi*(400/10000),2*%pi*(800/10000),0,0;

4

// Maximale Durchgangswelligkeit 0,5 dB

5

deltapass = (1-10**(-0.5/20)) ;

6

// minimale Sperrband-Dämpfung 50 dB

7

deltastop = 10**(-50/20);

8

// IIR-Filter berechnen, "lowpass", "butterworth":

9

[sos,gain,zeroes,poles] = eqiir('lp','butt',omega,deltapass,deltastop);

10

// Second order sections anzeigen:

11

sos

12

gain

13

//Frequenzgang berechnen und plotten:

14

zaehlerprodukt = prod(sos(2));

15

nennerprodukt = prod(sos(3));

16

frequenzgang = gain*abs(freq(zaehlerprodukt,nennerprodukt,exp(%i*(0:0.01:%pi))));

17

// frequenzgang = gain*abs(freq(zaehlerprodukt,nennerprodukt,exp(%i*(0:0.01:%pi))));

18

n = prod(size(frequenzgang));

19

//plot(20*log(frequenzgang(2:n))/log(10),(n*10000/(200*%pi)));

20

plot(20*log(frequenzgang (2:n))/log(10))

die 800 Hz Stopfrequenz und die dB-Werte habe ich erst mal willkürlich 
angenommen. Man müßte damit herumspielen, um auf ein Filter zweiter 
Ordnung zu kommen.

So ich habe das mal in Scilab ausprobiert, da fehlte eine eckige Klammer 
am Ende von omega=[...
Außerdem ergibt sich mit den Vorgaben ein IIR mit fast drei 
hintereinandergeschalteten second order sections. Hier ein IIR zweiter 
Ordnung:

1

// Lowpass-Filter für 400 Hz, Samplerate =10000Hz

2

// Eckfrequenzen 400Hz, 2000Hz:

3

omega=[2*%pi*(400/10000),2*%pi*(2000/10000),0,0];

4

// Maximale Durchgangswelligkeit 0,5 dB

5

deltapass = (1-10**(-0.5/20)) ;

6

// minimale Sperrband-Dämpfung 20 dB

7

deltastop = 10**(-20/20);

8

// IIR-Filter berechnen, "lowpass", "butterworth":

9

[sos,gain,zeroes,poles] = eqiir('lp','butt',omega,deltapass,deltastop);

10

// Second order sections anzeigen:

11

sos

12

gain

13

//Frequenzgang berechnen und plotten:

14

zaehlerprodukt = prod(sos(2));

15

nennerprodukt = prod(sos(3));

16

frequenzgang = gain*abs(freq(zaehlerprodukt,nennerprodukt,exp(%i*(0:0.01:%pi))));

17

// frequenzgang = gain*abs(freq(zaehlerprodukt,nennerprodukt,exp(%i*(0:0.01:%pi))));

18

n = prod(size(frequenzgang));

19

//plot(20*log(frequenzgang(2:n))/log(10),(n*10000/(200*%pi)));

20

plot(20*log(frequenzgang (2:n))/log(10))

21

22

23

Ergebnis:

24

25

sos  =

26

                      2          

27

            1 + 2z + z           

28

    --------------------------   

29

                              2  

30

    0.5275027 - 1.3735929z + z   

31

32

gain =  0.0384775

Die Durchgangsdämpfung ist allerdings noch unschön.

gegenüber der Formel hier:
http://www.mikrocontroller.net/articles/Digitalfilter_mit_ATmega#Von_der_.C3.9Cbertragungsfunktion_zum_Programm
sind oben die Buchstaben a und b vertauscht. Aus dem Ergebnis abgelesen 
folgt:

y(n) = 1*x(n) + 2*x(n-1) + 1*x(n-2) - 1.3735929*y(n-1) - 
0.5275027*y(n-2)

Hi Christoph,

vielen Dank, dass du dir die Mühe gemacht hast es zu berechnen,
ich bin auch gerade dabei es zu verstehen.

So langsam aber sicher, fange ich an, dort durchzublicken!

Was meinst du denn, mit der Durchgangsdämpfung ?

"gain =  0.0384775" das ist die Durchgangsverstärkung, genauer eine 
Dämpfung. Ein Sinus (< 400 Hz) der Amplitude "1" am Eingang kommt nur 
mit der Amplitude 0,038 wieder heraus. Das stört vor allem, wenn man nur 
mit wenigen Bit Auflösung arbeitet.

Die besten Filtereigenschaften erreicht man ungefähr mit 
Halbbandfiltern, das heißt, die Filterflanke ist symmetrisch zur halben 
Nyquistfrequenz oder einem Viertel der Abtastrate. Das wären hier 2500 
Hz. Das IIR-Filter hat ja trotz des Namens nur ein sehr begrenztes 
"Erinnerungsvermögen", sonst hat man einen Oszillator gebaut. Je 
niedriger die Frequenz im Verhältnis zur Abtastrate, desto schlechter 
werden die Eigenschaften.

Zweiter Ordnung und noch eine "gemütliche" Filtercharaktristik wie 
Bessel oder Butterworth führt leider zu dem langsamen Abfall der Kurve. 
Bei der fünffachen Grenzfrequenz 2000Hz sind es hier erst 20 dB 
Sperrdämpfung.
Wenn die größeren Phasenänderungen nahe der Grenzfrequenz nicht stören, 
würde ich erst mal ein "elliptic" probieren, das hat die größte 
Steilheit bei gleicher Ordnungszahl.