Hallo :)
Ich versuche grade einen Digitalen Bandpass in C zu implementieren.
Leider klappt das nicht ganz so wie gewünscht.
Ich versuche folgende Schaltung nachzubauen:
1 | Uin ____ Uout
|
2 | o------|____|---+----+--o
|
3 | R | |
|
4 | L _|_C
|
5 | L ___
|
6 | L |
|
7 | | |
|
8 | +----+
|
9 | |
|
10 | _|_
|
11 | -
|
Die Übertragungsfunktion für das System ist:
wobei meine Resonanzfrequenz
ist und meine Bandbreite
Damit habe ich eine Null bei
und 2 Pole bei
Jetzt will ich diesen Filter in die Z-Ebene bringen.
Dieser Seite [1] habe ich die Transformation
entnommen, wobei T die Periodendauer meiner Abtastung ist.
Damit bekomme ich für die Nullstelle
und für die erste Polstelle:
Die zweite Polstelle ist entsprechend komplex konjugiert.
Nun multipliziere ich den Nenner aus:
Jetzt erhalte ich für mein ganzes Polynom:
Das nurnoch kürzen:
Erstmal, ist das soweit richtig?
Als nächstes habe ich einen konkreten Fall gewählt, und zwar
und
und das ganze in C++ implementiert. Als Umsetzung der Filterstruktur
habe ich die Direktform 2 gewählt [2].
Hier der dazugehörige Code:
1 | #define _USE_MATH_DEFINES
|
2 | #include <math.h>
|
3 | #include <iostream>
|
4 |
|
5 | int main()
|
6 | {
|
7 | double w0 = 100 * 2 * M_PI; // 100Hz
|
8 | double w_band = 2 * M_PI * 10; // 10Hz -3db
|
9 | double samplerate = 1000; // abtastrate in Hz
|
10 |
|
11 | double w_res_sq = w0 * w0;
|
12 |
|
13 | double pol_real = -w_band / 2; // realteil des pols
|
14 | double pol_img = sqrt(w_res_sq - w_band * w_band / 4); // imaginärteil des pols
|
15 | double exp_real = pol_real / samplerate; // mit T skaliert
|
16 | double exp_img = pol_img / samplerate;
|
17 | // e^(a+bi) = e^a * e^(bi) = e^a * (cos(b)+i*sin(b))
|
18 | double x = exp(exp_real) * cos(exp_img); // lage des pols in der z-ebene
|
19 | double y = exp(exp_real) * sin(exp_img);
|
20 |
|
21 | // zähler
|
22 | double b0 = 0;
|
23 | double b1 = 1;
|
24 | double b2 = -1;
|
25 | // nenner
|
26 | double a0 = 1;
|
27 | double a1 = -2 * x;
|
28 | double a2 = x * x + y * y;
|
29 |
|
30 | double val0 = 0, val1 = 0, val2 = 0;
|
31 | double time = 0;
|
32 | for (unsigned i = 0; i < 1000; ++i)
|
33 | {
|
34 | double newval = sin(time * w0);
|
35 | time += 1 / samplerate;
|
36 | val2 = val1;
|
37 | val1 = val0;
|
38 | val0 = newval + b1 * val1 + b2 * val2;
|
39 |
|
40 | double out = a0 * val0 + a1 * val1 + a2 * val2;
|
41 |
|
42 | std::cout << "in = " << newval << std::endl;
|
43 | std::cout << "out = " << out << std::endl;
|
44 | }
|
45 | }
|
Die Idee war dass ich ein Sinussignal auf den Eingang gebe, welches
exakt die Resonanzfrequenz des Filters hat. An sich sollte der Filter
langsam anfangen zu schwingen und am Schluss identisch zum Signal
aussehen. Leider nimmt der Filter aber werte größer als 1 an, was bei
diesem Eingangssignal nicht möglich ist. Daher scheine ich wohl irgentwo
ein Fehler zu haben :(
Die Filterkoeffizienten sind:
1 | b0 = 0
|
2 | b1 = 1
|
3 | b2 = -1
|
4 | a0 = 1
|
5 | a1 = -1.5688869346931174
|
6 | a2 = 0.93910136742429273
|
Grüße,
Resonanzkatastrophe
Quellen:
[1]
http://www.dummies.com/how-to/content/working-across-domains-example-take-the-rc-lowpass.html
[2]
http://de.wikipedia.org/wiki/Filter_mit_unendlicher_Impulsantwort#Direkt-Form_2_.28DF2.29
