Forum: Digitale Signalverarbeitung / DSP / Machine Learning Filterkoeffizienten für IIR Butterworth 3. Ordung berechnen?


von Jürgen H. (t-5)


Lesenswert?

Servus!

Ich habe einen Satz LADSPA-Plugins für eine 
Lautsprecher-Frequenzweichen-Software unter Linux geschrieben:

https://t-5.eu/hp/Software/ladspa-t5-plugins/

https://t-5.eu/hp/Software/Pulseaudio%20Crossover%20Rack/

Auf Anfrage würde ich nun gerne noch Butterworth Hoch- und Tiefpässe 3. 
Ordnung hinzufügen und bräuchte dafür Berechnungsformeln für die 
Filterkoeffizienten wie z.B. für die 2. Ordnung:
1
    def _coeffsBA(f, samplerate):
2
        """ return 12db/oct high pass biquad coeffs for given frequency """
3
        w0 = 2 * pi * f / samplerate
4
        alpha = sin(w0) / 2 / 0.7071067811865476 # Butterworth characteristic, Q = 0.707...
5
        cs = cos(w0)
6
        norm = 1 / (1 + alpha)
7
        b0 = (1 + cs) / 2 * norm
8
        b1 = -1.0 * (1 + cs) * norm
9
        b2 = b0
10
        a1 = -2 * cs * norm
11
        a2 = (1 - alpha) * norm
12
        return (b0, b1, b2), (1, a1, a2)

Für erster bzw. zweiter Ordnung war einiges zu finden. Für dritte 
Ordnung leider gar nix.

Hat da jemand was für mich?

VG Jürgen

von Dergute W. (derguteweka)


Lesenswert?

Moin,

Kannst dir ja mal das hier angucken:

https://octave.sourceforge.io/signal/

in dem Package gibts eine funktion namens butter; die macht, das was du 
gerne haettest.
Wird nur bissl aufwendiger, weil ein Filter mit Ordnung >2 
ueblicherweise in Teilfilter 2. Ordnung und ggf. noch eines 1. Ordnng 
aufgeteilt werden. Das sind dann noch extra Funktionen.

Gruss
WK

von Jürgen H. (t-5)


Lesenswert?

Dergute W. schrieb:
> Kannst dir ja mal das hier angucken:
>
> https://octave.sourceforge.io/signal/
>
> in dem Package gibts eine funktion namens butter; die macht, das was du
> gerne haettest.

Danke! Ich werd mir das morgen mal zu Gemüte führen. Ich hab auch schon 
in die Source des python-Pakets "scipy.signal" geschaut, wo es auch eine 
funktion butter() gibt. Leider etwas aufwändig zu reverse-engineeren, 
weil da Zwischenformen verwendet werden, die mir momentan gar nix sagen 
(Stichwort zpk aka pole-zero). Jetzt gehts dann erstmal ins Bett, morgen 
dann Code-Studie bei octave :)

> Wird nur bissl aufwendiger, weil ein Filter mit Ordnung >2
> ueblicherweise in Teilfilter 2. Ordnung und ggf. noch eines 1. Ordnng
> aufgeteilt werden. Das sind dann noch extra Funktionen.

Kannst Du mir nen Hinweis geben, wo ich was darüber lesen kann? Formeln 
für die Filter erster und zweiter Ordnung hab ich ja schon zur Hand. 
Aber ich gehe mal stark davon aus, dass da mit Verschiebung von 
Filterfrequenzen und -güte gearbeitet wird, wenn die Filter kaskadiert 
werden. Ansonsten kommt ja nicht eine Dämpfung von 3dB bei der 
Grenzfrequenz raus sondern 6dB, oder?

Alles in allem sehr sehr interessant das Thema. Ich lerne fast täglich 
:) Schließlich hab ich das Projekt erst vor ein paar Wochen angefangen 
mit praktisch Null Wissen bzgl C-Programmierung und drüben auf 
diyaudio.com setzt das Programm wohl schon jemand aktiv für seine 
Lautsprecher ein...

VG Jürgen

von Dergute W. (derguteweka)


Lesenswert?

Moin,

Fuer den Anfang willst du da wohl gleich ganz schoen dicke Bretter 
bohren :-)
Ist jetzt dummerweise so, dass die "Formeln" fuer Filter hoeherer 
Ordnung nicht mehr ganz so simpel sind, wie die fuer 1 und 2. Ordnung.
Also sowas wie Pol/Nullstellen usw. gehoeren da zum Handwerkszeug.
Vielleicht hilft das:

https://www.mikroe.com/ebooks/digital-filter-design/introduction-iir-filter

Oder wahrscheinlich jedes andere gute Buch ueber analoges Filterdesign 
gefolgt von einem anderen guten Buch ueber Digitalfilter...

Gruss
WK

von Jürgen H. (t-5)


Lesenswert?

Dergute W. schrieb:
> Moin,
>
> Fuer den Anfang willst du da wohl gleich ganz schoen dicke Bretter
> bohren :-)
> Ist jetzt dummerweise so, dass die "Formeln" fuer Filter hoeherer
> Ordnung nicht mehr ganz so simpel sind, wie die fuer 1 und 2. Ordnung.
> Also sowas wie Pol/Nullstellen usw. gehoeren da zum Handwerkszeug.

Moinsen!

für das erste Projekt, den Pulseaudio Parametric Equalizer, hab ich ja 
noch quasi vorgebohrte Bretter verwenden können, nämlich aus dem Audio 
EQ Cookbook. Für die Crossover-Filter wurde es schon etwas schwieriger, 
komischerweise war 2. Ordnung einfacher zu finden als 1. Ordnung, die 
ich mir etwas umständlicher herleiten musste. Und nun, ja da muss ich 
mir wohl eine Bohrmaschine anschaffen, um im Bild zu bleiben :)))

> Vielleicht hilft das:
>
> https://www.mikroe.com/ebooks/digital-filter-design/introduction-iir-filter

Danke, ich werde lesen und berichten!

VG Jürgen

von Gerb (Gast)


Lesenswert?

Jürgen H. schrieb:
> Kannst Du mir nen Hinweis geben, wo ich was darüber lesen kann? Formeln
> für die Filter erster und zweiter Ordnung hab ich ja schon zur Hand.

Ich bin für meine Filterberechnungen sehr gut mit dem Tietze Schenk 
gefahren (Halbleiter-Schaltungstechnik: erst Kapitel über analoge, dann 
das über digitale Filter lesen) mit Tabellen, um Koeffizienten zu 
berechnen.

Gerhard

von Daniel (Gast)


Lesenswert?

Der Sourcecode von dfcgen-gtk sieht sehr aufgeräumt aus.
Es wird aber für einige Sachen auf die GSL zurückgegriffen.

von Edi M. (Gast)


Lesenswert?

Dergute W. schrieb:
> Wird nur bissl aufwendiger, weil ein Filter mit Ordnung >2
> ueblicherweise in Teilfilter 2. Ordnung und ggf. noch eines 1. Ordnng
> aufgeteilt werden. Das sind dann noch extra Funktionen.

Ist das denn dann überhaupt noch analytisch beschreib- und lösbar oder 
müsste man eine Iteration vornehmen, um die Lösung zu finden?

von Dergute W. (derguteweka)


Lesenswert?

Moin,

Edi M. schrieb:
> Ist das denn dann überhaupt noch analytisch beschreib- und lösbar oder
> müsste man eine Iteration vornehmen, um die Lösung zu finden?

Ich kenn' kein analytisches Verfahren, um Nullstellen in Polynomen >2. 
Ordnung zu finden. In der Schulmathematik wird iirc dann immer eine 
"simple" Nullstelle geraten. Das haut natuerlich im real-life nicht so 
hin.


Gruss
WK

von Detlef _. (detlef_a)


Lesenswert?

>>Ich kenn' kein analytisches Verfahren, um Nullstellen in Polynomen >2.
>>Ordnung zu finden.

Für 3. und 4. Ordnung gibts die noch, darüber nicht mehr:
https://de.wikipedia.org/wiki/Cardanische_Formeln

Numerische Berechnung in Matlab heißt 'roots'.

math rulez!
Cheers
Detlef

von Dergute W. (derguteweka)


Lesenswert?

Moin,

Sapperlott, den Cardano hab' ich noch nicht gekannt. Wieder was gelernt!

Gruss
WK

von Helmut S. (helmuts)


Lesenswert?

Für das Butterworth-Filter muss man gar keine Gleichung lösen.
https://de.wikipedia.org/wiki/Butterworth-Filter
Dafür reicht ein Rechner der cos() kann.

3. Grad

F(p) = 1/((1+p)*(1+p+p^2))

F(p) = 1/(1+2*p+2*p^2+p^3)

: Bearbeitet durch User
von Edi M. (Gast)


Lesenswert?

Hey, ein paar neue Denkanstösse. Sehr gut. Man sieht wieder mal, dass 
man nicht genug Mathe gehabt haben kann!

von Jürgen H. (t-5)


Lesenswert?

Helmut S. schrieb:
> Für das Butterworth-Filter muss man gar keine Gleichung lösen.
> https://de.wikipedia.org/wiki/Butterworth-Filter
> Dafür reicht ein Rechner der cos() kann.
>
> 3. Grad
>
> F(p) = 1/((1+p)*(1+p+p^2))
>
> F(p) = 1/(1+2*p+2*p^2+p^3)

Ich hab's immer noch ned ganz kapiert, was da genau gemacht wird, aber 
nach lesen der scipy.signal-sourcen, von 
https://www.mikroe.com/ebooks/digital-filter-design/introduction-iir-filter 
und https://www.dsprelated.com/showarticle/1119.php bin ich schon 
ewtas schlauer.

Wenn ich o.g. kubische Gleichung lösen kann, bekomme ich die 3 Pole der 
Übertragungsfunktion als komplexe Zahlen, ist das soweit korrekt? Und 
das für den analogen Filter. Auch korrekt soweit?

Was mir also noch fehlt ist das Prewarping der Grenzfrequenz im 
diskreten Filter, das ist (ausnahmsweise) easy.

Und dann die Bilineare Transformation in die Zeitdiskrete Form. bei 
letzterem hakt's in C zumindest noch gewaltig. Ich versuch' mich grad an 
der GSL, aber das wird noch etwas dauern, bis ich das alles kapiere. In 
der Zwischenzeit wäre ich schon mal dankbar für Input zu meinen Zeilen 
davor...

von Jürgen H. (t-5)


Lesenswert?

Keiner mehr was zu sagen? :)

von Jürgen H. (t-5)


Lesenswert?

Falls es jemanden interessiert, ich hab den Code für die Funktion 
butter() in scipy.signal reverse engineered und vereinfacht für den Fall 
N=3 (3. Ordnung) und dann soweit möglich in C umgesetzt.

Vorläufiges Resultat: https://pastebin.com/sPVtgB1s

Was ich noch brauche ist der Code von numpy.poly(), das aus den 
gegebenen Nullstellen die Faktoren eines Polynoms berechnet (Die Stelle 
mit "NEED HELP" im oben gezeigten Code.

Kann hier evtl. jemand helfen?

von Detlef _. (detlef_a)


Lesenswert?

Ich verstehe nicht, was Du tust. In dem von Dir zitierten dsprelated 
Artikel steht doch alles was Du willst als Kochrezept, da gibts' nix was 
man reengineeren müßte.

Ausmultiplizieren von Nullstellen geht über eine Faltung, bei Matlab 
conv:

conv Convolution and polynomial multiplication.
    C = conv(A, B) convolves vectors A and B.  The resulting vector is
    length MAX([LENGTH(A)+LENGTH(B)-1,LENGTH(A),LENGTH(B)]). If A and B 
are
    vectors of polynomial coefficients, convolving them is equivalent to
    multiplying the two polynomials.

Das machste einfach für alle Nullstellen. Kann man sich aber auch 
schlicht selber zurechtlegen, wenn man sich vorstellt, wie man an ein 
vorhandenes Polynom ein weiteres Polynom 1. Ordnung mit einer Nullstelle 
ranmultipliziert.

math rulez
Cheers
Detlef

von Jürgen H. (t-5)


Lesenswert?

Jo, ich hab nun quasi den matlab bzw scipy-code in C implementiert.

Wer's sehen will:
https://gitlab.com/t-5/ladspa-t5-plugins/blob/master/src/coeffs.h

Scheint soweit alles zu funktionieren, wenn's auch vermutlich nicht grad 
elegant gelöst is mit meinem C-Programmierstil :)

Danke an alle, die geholfen haben!

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
Noch kein Account? Hier anmelden.