Design eines Bandpasses

[Resonanzkatastrophe (Gast)]

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:

 Uin     ____          Uout
 o------|____|---+----+--o
           R     |    |
                 L   _|_C
                 L   ___
                 L    |
                 |    |
                 +----+
                      |
                     _|_
                      -

Die Übertragungsfunktion für das System ist:

$$G(s)=\frac{s/RC}{s^2+s/RC+1/LC}$$

wobei meine Resonanzfrequenz

$$w_{res}=1/sqrt(LC)$$

ist und meine Bandbreite

$$w_{band}=1/RC$$

Damit habe ich eine Null bei

$$s=0$$

und 2 Pole bei

$$s=\sqrt{\frac{w_{band}^2}{4}-w_{res}^2}-\frac{w_{band}}{2}$$

Jetzt will ich diesen Filter in die Z-Ebene bringen.

Dieser Seite [1] habe ich die Transformation

$$z=e^{sT}$$

entnommen, wobei T die Periodendauer meiner Abtastung ist.

Damit bekomme ich für die Nullstelle

$$z=1$$

und für die erste Polstelle:

$$s_{real}=-\frac{w_{band}}{2}$$

$$s_{img}=\sqrt{w_{res}^2-\frac{w_{band}^2}{4}}$$

$$z_{real}=e^{s_{real}T} \cdot \cos(s_{img} \cdot T)$$

$$z_{img}=e^{s_{real}T} \cdot \sin(s_{img} \cdot T)$$

Die zweite Polstelle ist entsprechend komplex konjugiert.

Nun multipliziere ich den Nenner aus:

$$(z-(z_{real}+j z_{img}))(z-(z_{real}-j z_{img})) = z^2 - 2 z_{real} z + (z_{real}^2 + z_{img}^2)$$

Jetzt erhalte ich für mein ganzes Polynom:

$$H(z)=\frac{z-1}{z^2 - 2 z_{real} z + (z_{real}^2 + z_{img}^2)}$$

Das nurnoch kürzen:

$$H(z)=\frac{z^{-1}-z^{-2}}{1 - 2 z_{real} z^{-1} + (z_{real}^2 + z_{img}^2) z^{-2}}$$

Erstmal, ist das soweit richtig?

Als nächstes habe ich einen konkreten Fall gewählt, und zwar

$$f_{res}=100Hz$$

und

$$f_{band}=10Hz$$

und das ganze in C++ implementiert. Als Umsetzung der Filterstruktur
habe ich die Direktform 2 gewählt [2].

Hier der dazugehörige Code:

#define _USE_MATH_DEFINES
#include <math.h>
#include <iostream>

int main()
{
  double w0 = 100 * 2 * M_PI; // 100Hz
  double w_band = 2 * M_PI * 10; // 10Hz -3db
  double samplerate = 1000; // abtastrate in Hz

  double w_res_sq = w0 * w0;
  
  double pol_real = -w_band / 2; // realteil des pols
  double pol_img = sqrt(w_res_sq - w_band * w_band / 4); // imaginärteil des pols
  double exp_real = pol_real / samplerate; // mit T skaliert
  double exp_img = pol_img / samplerate;
  // e^(a+bi) = e^a * e^(bi) = e^a * (cos(b)+i*sin(b))
  double x = exp(exp_real) * cos(exp_img); // lage des pols in der z-ebene
  double y = exp(exp_real) * sin(exp_img);

  // zähler
  double b0 = 0;
  double b1 = 1;
  double b2 = -1;
  // nenner
  double a0 = 1;
  double a1 = -2 * x;
  double a2 = x * x + y * y;

  double val0 = 0, val1 = 0, val2 = 0;
  double time = 0;
  for (unsigned i = 0; i < 1000; ++i)
  {
    double newval = sin(time * w0);
    time += 1 / samplerate;
    val2 = val1;
    val1 = val0;
    val0 = newval + b1 * val1 + b2 * val2;

    double out = a0 * val0 + a1 * val1 + a2 * val2;

    std::cout << "in = " << newval << std::endl;
    std::cout << "out = " << out << std::endl;
  }
}

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:

b0 = 0
b1 = 1
b2 = -1
a0 = 1
a1 = -1.5688869346931174
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

[abcd (Gast)]

Ich habe deine Übertragungsfunktion nur mal kurz durch Matlab gejagt. 
Deine Filterkoeffizienten stimmen auf jeden Fall, der Frequenzgang hat 
auch die gewollte Spitze bei 100Hz und zeigt das von dir erwartete 
Verhalten.

Das Problem ist, dass du deinen Filter noch auf eine bestimmte 
Verstärkung normieren musst. Das wird in deinem Link [1] auch unter 
"Fine-tune the design for implementation" erklärt.

In deinem Fall wäre es wohl am Einfachsten sich die Verstärkung im 
kontinuierlichen Filter bei der Resonanzfrequenz und im diskreten Filter 
bei der Resonanzfrequenz anzuschauen und dann den Gainfactor 
entsprechend zu wählen, dass es passt.

[Resonanzkatastrophe (Gast)]

Hallo abcd

Du hast Recht, der Verstärkungsfaktor fehlt tatsächlich. Aber auch nach 
einfügen selbigens hatte der Filter ein komisches Verhalten, die 
eigentliche Ursache war im Code. Die Filterkoeffizienten waren 
vertauscht.

vorher:

val2 = val1;
val1 = val0;
val0 = newval + b1 * val1 + b2 * val2;

double out = a0 * val0 + a1 * val1 + a2 * val2;

nachher:

val2 = val1;
val1 = val0;
val0 = newval - a1 * val1 - a2 * val2;

double out = (b0 * val0 + b1 * val1 + b2 * val2) / scale;

Zum einen waren nenner und zählerkoeffizienten falsch und zum anderen 
muss ich subtrahieren statt addieren. So funktioniert der Filter wie 
gewünscht, nur weiß ich leider nicht warum das sorum sein muss. Bei der 
Abbildung auf Wikipedia waren die koeffizienten schließlich anders rum.

Viele Grüße,
Resonanzkatastrophe

[Freddy (Gast)]

Resonanzkatastrophe schrieb:
> Bei der
> Abbildung auf Wikipedia waren die koeffizienten schließlich anders rum

In der Literatur ist man sich uneins, ob b die Zähler- und a die 
Nenner-Koeffizienten oder andersrum sind.

Bei dem Wikipedia-Artikel haben sie bei den Gleichungen b im Zähler so 
wie du, aber mit den Bildern waren sie dann ein bisschen inkonsequent:

$$\alpha \hat{=} b$$

$$\beta \hat{=} \mathbf{-}a$$

Ersteres erklärt dein Vertauschen der Koeffizienten, zweiteres das 
nötige Subtrahieren.