Implementierung digitaler Filter

[Bert (Gast)]

Hallo zusammen,

befasse mich gerade ein wenig mit digitaler Signalverarbeitung.

Ich habe gelesen, dass für Audioverarbeitung oft FIR-Filter in Lattice 
Struktur eingesetzt werden. Meine Frage lautet nun: was hat die 
Lattice-Struktur gegenüber einem 'normalen' FIR-Filter für einen 
Vorteil? Warum eignet sie sich angeblich besonders gut für Audio? Sind 
die Reflexionskoeffizienten irgendwie weniger anfällig auf 
Quantisierungsfehler (Runden aufgrund von Fixpoint-Implementierung)?

Wenn ich eine Audiosignalverarbeitung implementieren möchte und dazu ein 
solches Lattice-Filter verwende, wie implementiere ich das dann 
SW-Mässig auf meinem DSP/uC ? Oder würde man zur 
Software-Implementierung einfach ein gewöhnliches FIR nehmen?

Dieselbe Frage stelle ich mir auch bei IIR-Filtern: dort wird immer 
gesagt, dass man die üblicherweise in eine SOS-Form und einen Gain 
zerlegt. Wird bei der Software-Implementierung tatsächlich dann auch mit 
einzelnen SOS-Filtern gerechnet und am Schluss kommt noch der Gain 
hinzu, oder wie soll ich mir das vorstellen?

Den Sinn dieser verschiedenen Filtertopologien sehe ich schon ein, wenn 
man das auf einem FPGA realisieren würde. Aber wenn ich es in SW auf 
einem DSP mache habe ich fast das Gefühl, der Aufwand, 5 SOS-Filter 
nacheinander zu berechnen, ist grösser, als wenn ich ein einziges Filter 
10. Ordnung rechne. Und dazu habe ich dann auch mehr Rundungsprobleme. 
Beim Lattice denke ich mir das ähnlich - da ist doch ein simples FIR 
einfacher, als wenn ich mir die Mühe mache, in SW einen solchen Lattice 
zu beschreiben! oder täusche ich mich da?

Gruss!

[Johannes E. (cpt_nemo)]

Bert schrieb:
> Aber wenn ich es in SW auf
> einem DSP mache habe ich fast das Gefühl, der Aufwand, 5 SOS-Filter
> nacheinander zu berechnen, ist grösser, als wenn ich ein einziges Filter
> 10. Ordnung rechne. Und dazu habe ich dann auch mehr Rundungsprobleme.

Bei einem Filter 10. Ordnung müssen die Filterkoeffizienten relativ 
genau eingehalten werden, damit das Filter die gewünschte Eigenschaften 
hat.
Mit 5 einzelnen Filtern 2. Ordnung hat man eine größere Toleranz 
gegenüber Rundungsfehlern bei der Berechnung bzw. Quantisierung der 
Filterkoeffizienten.

Speziell bei IIR-Filtern hoher Ordnung können durch die Rückkopplung 
auch sehr kleine Fehler eine große Wirkung haben.

[Bert (Gast)]

Hallo Johannes

ja, das Problem ist mir bekannt. Darum teilt man ein grösseres Filter ja 
auch in SOS-Teilfilter auf. Sind die Teilfilter für sich stabil, dann 
ist es auch das aus der Serieschaltung resultierende Gesamtfilter. Und 
für ein Filter 2. Ordnung lässt es sich noch verhältnismässig einfach 
und schnell feststellen, ob es stabil ist.

Aber wie gesagt, wie realisiere ich dies in C auf einem DSP / 
Mikrocontroller? Dass man in einem FPGA einfach 5 Filter 2. Ordnung 
baue, um mein IIR-Filter 10.Ordnung zu erhalten, ist natürlich klar, und 
ausser ein paar zusätzlichen LEs oder Slices bringt das keinen 
Zusatzaufwand. Aber wie sieht es auf dem Mikrocontroller aus? Dort kann 
ich ja nicht die 5 Filter parallel rechnen.

Macht man es etwa so (Pseudocode):

y1[k] = b01*u[k] + b11*u[k-1] + b12*u[k-2] - a11*y1[k-1] - a21*y1[k-2];

y2[k] = b02*y1[k] + b12*y1[k-1] + b22*y1[k-2] - a12*y2[k-1] - a22*y2[k-2];
...

man rechnet also einfach mit ein "paar" Zwischenresultaten? ist das so 
nicht super-duper ineffizient?

[Johannes E. (cpt_nemo)]

Bert schrieb:
> man rechnet also einfach mit ein "paar" Zwischenresultaten? ist das so
> nicht super-duper ineffizient?

Es ist schon etwas mehr Rechenaufwand, aber so extrem ist das dann auch 
nicht:

Ein Filter 2. Ordnung hat 5 Multiplikationen, wenn man 5 solcher Filter 
durchrechnet, sind das 25 Multiplikationen.

Ein Filter 10. Ordnung hat 21 Multiplikationen, also 4 weniger als mit 
den einzelnen Filtern.

Bei der Anzahl der Additionen ist es ähnlich. Die Frage bei so etwas ist 
nicht primär die Effizienz, sondern erst mal ist wichtig, dass es 
überhaupt funktioniert. Wenn ein Filter 10. Ordnung etwas effizienter 
ist als 5 Filter 2. Ordnung, aber durch die Quantisierungsfehler 
instabvil ist, hat man ja auch nichts gewonnen.

[Bert (Gast)]

Hallo Johannes,
ah okay, man macht es also wirklich genauso, wie man auch das 
Blockschaltbild zeichnet - indem effektiv mehrere Filter 2. Ordnung 
hintereinander kommen. Danke für die Auskunft.

Hast du soetwas schon einmal so realisiert?


Zum Lattice-Filter kannst du auch etwas sagen?

Gruss & Dank

[Bert (Gast)]

Hallo,

Nachtrag zu meinem vorherigen Post: wäre folgendes eine adäquate 
Implementierung? u ist der Puffer für die alten Eingangswerte. Der ist 
einfach. Dann, y ist ein zweidimensionales Array, dessen erste Dimension 
die Sektion angibt, und die zweite Dimension sind dann einfach die alten 
Werte.

static float sos_iir(float in)
{
  static float u[3] = { 0 };
  static float y[5][3] = { 0 };
  u[0] = in;
  
  
  
  u[2] = u[1]; u[1] = u[0];
  y[0][2] = y[0][1]; y[0][1] = y[0][0];
  y[0][0] = b[0][0]*u[0] +
            b[0][1]*u[1] +
            b[0][2]*u[2] -
            a[0][1]*y[0][1] -
            a[0][2]*y[0][2];
  
  for(i = 1; i < 5; i++)
  {
    y[i][2] = y[i][1]; y[i][1] = y[i][0];
    y[i][0] = b[i][0]*y[i-1][0] +
              b[i][1]*y[i-1][1] +
              b[i][2]*y[i-1][2] -
              a[i][1]*y[i][1] -
              a[i][2]*y[i][2];
  }
  
  return y[4][0];
}

[Johannes E. (cpt_nemo)]

Bert schrieb:
> Nachtrag zu meinem vorherigen Post: wäre folgendes eine adäquate
> Implementierung?

Ich denke schon, dass es so funktioniert.

Es gibt bei so etwas natürlich immer unterschiedliche Lösungswege. Die 
Frage ist, welche Eigenschaften dir wichtig sind bzw. auf was man den 
Code optimieren möchte (Code-Größe, RAM-Verbrauch, Brechnungsdauer).

Man kann es z.B. auch so machen, dass man eine Funktion schreibt, die 
einen Filter 2. Ordnung durchrechnet, die Speicher-Arrays werden als 
Funktions-Prameter übergeben.
Diese Funktion kann man dann 5 mal hintereinander aufrufen und übergibt 
jeweils den passenden Pointer. Dadurch sollte die Code-Größe kleiner 
werden.

Wenn man die Filter in der "Direkt-Form 2" 
(http://de.wikipedia.org/wiki/Filter_mit_unendlicher_Impulsantwort#Direkt-Form_2_.28DF2.29) 
implementiert, wird weniger RAM benötigt.

Ein wichtiger Punkt bei digitalen Filtern ist das Zahlenformat. Du hast 
in deinem Beispiel Fließkomma-Zahlen verwendet; das ist auf einem 
Controller ohne Fließkomma-Einheit sehr ungünstig; da verwenet man 
normalerweise Festkomma-Zahlen.

Auf einem DSP gibt es meistens spezielle MAC-Befehle (Multiply and 
accuulate), die man aber nicht so einfach aus C-Code heraus aufrufen 
kann. Hier muss man dann entweder den Filter in Assembler programmieren 
oder eine für diese CPU optimierte Bibliothek verwenden.


Du siehst, es ist also nicht so einfach, "den besten" Algotithmus 
auszuwählen, das hängt von sehr vielen Randbedinungen ab.

[Bert (Gast)]

Hallo,

vielen Dank für deine Ausführungen. Das Filter ansich soll später auf 
einem MSP430 laufen. Aber da ich dazu die Hardware noch nicht habe, 
werde ich zunächst das Cortex M4 Discoveryboard benutzen. Der dort 
verbaute Microcontroller hat auch eine FPU, sodass ich getrost 
Fliesskomma verwenden kann. Natürlich muss ich es dann für den MSP in 
ein gescheites Zahlenformat ummünzen, das war mir klar.

Noch eine Frage:
In Matlab verwende ich den Befehl tf2sos, um mein Filter zu bestimmen:

http://www.mathworks.ch/ch/help/signal/ref/tf2sos.html

wie der Algorithmus funktioniert, ist mir klar, und wenn ich es von Hand 
durchrechne, komme ich auch auf das selbe Ergebnis wie Matlab. Das 
scheint also zu funktionieren. SOS ist eine Matrix mit 6 Spalten (b0 b1 
b2 1 a1 a2) und 5 Zeilen (bei dem Filter 10. Ordnung). G ist der Gain, 
welcher natürlich Skalar ist.

Wenn ich in Matlab allerdings mti dem FDATOOL ein Filter baue, dann 
bekomme ich als Resultat die SOS Matrix und einen Gain-*Vektor* - der 
Gain wird da auf die einzelnen Sektionen aufgeteilt, aber nicht jede 
Sektion hat gleich viel Gain. Was hat das für einen Sinn? Muss ich das 
auch machen oder kann ich einfach einen Gesamtgain verwenden?

[Ingenieur (Gast)]

Die Gainproblematik resultiert aus den unterschiedlichen Gains, der 
filter-Koeffs, man spart sich da einfach Zwischenskalierungen auf 1. 
Einen anderen Grund sehe ich da nicht.

Bei Audio musst Du bedenken, dass dort oft mit 16Bit und mehr an 
Genaugkeit gearbeitet werden soll. Ein Filter auf 90dB glatt zu bekommen 
ist aufwändig.

Anders herum sind Audio-Apss oft billig realisiert, der Kosten wegen. 
Ich würde also auch davon ausgehen, dass die Filter, die Du dort 
vorfindest, eine Kompromiss sind.

[Johannes E. (cpt_nemo)]

Bert schrieb:
> Muss ich das
> auch machen oder kann ich einfach einen Gesamtgain verwenden?

Bei Fließkomma-Zahlen ist es vermutlich egal, wie du es machst.
Bei Festkomma-Zahlen sollte man drauf achten, dass die Zahlenbereiche 
möglichst optimal ausgenutzt werden, da kann es schon sinnvoll sein, den 
Gain geschickt aufzuteilen.

[Bert (Gast)]

Okay. Also dann wäre es bei Fixpoint ggf. besser, den Gain aufzuteilen. 
Ich würde das dann so machen:

Für jedes Teilfilter wird der DC-Gain ausgerechnet. Und den Gesamtgain 
teile ich dann so auf jedes Teilfilter auf, dass bei jedem Teilfilter

(DC Gain) * (aufgeteilter Gain)

maximal wird, ja? würdest du auch so machen?

Ich habe übrigens jetzt gerade meine Implementierung des genannten 
Filters 10. Ordnung in C implementiert. Und mal laufen lassen. Es 
handelt sich um den Code, den ich weiter oben gepostet habe. Man kann 
direkt die SOS-Matrix von Matlab in meinen Code rein kopieren, die 
Filterlänge angeben, und fertig, das ergebnis (Impuls- und 
Sprungantwort) stimmt exakt mit Matlab überein. (Jetzt ist es in float 
implementiert, der Cortex M4 dürfte damit gar kein Problem haben.)
Den Code kann vielleicht sonst noch jemand brauchen, darum habe ich ihn 
angehangen!

Jetzt befasse ich mich noch mit dem Lattice Filter. Ich habe dazu mal 
etwas geschrieben, allerdings habe ich bei dem Code einen Denkfehler 
gemacht, er funktioniert nicht. Ich meine, es sollte beim Lattice Filter 
möglich sein, den Filter Output zu berechnen, ohne dass ich die Samples 
rumkopieren muss. Beim FIR muss ich ja immer etwas in der Art

u[5] = u[4]; u[4] = u[3]; u[3] = u[2]

usw. schreiben; ich würde mal meinen, beim Lattice liesse sich das 
umgehen, was natürlich  interessant wäre. Aber so oder so, bekomme ich 
es nicht gebacken :-)
Vorerst würde ich mich mit einem Lattice FIR begnügen:
http://de.wikipedia.org/wiki/Latticefilter
Vielleicht kommt ja noch Karl Heinz Buchegger hier vorbei und kann mir 
verraten, wo mein Denkfehler liegt :-)

#include <stdio.h>


#define L 6
static const float K[L] =
{

   100.0000e-003, 464.6465e-003, 846.8445e-003, -972.9879e-003, 23.0590e+000, 965.1336e-003,
};

/* delay elements */
static float f[L];
static float g[L];


static float u_old = 0.0f;

static float lattice_fir(float u)
{
  int i;
  
  /* special case: 1st stage */
  f[0] = u + K[0]*u_old;
  g[0] = u_old + K[0]*u;
  
  u_old = u;
  
  /* handle all other stages */
  for(i = 1; i < L; i++)
  {
    f[i] = f[i-1] + K[i]*g[i-1];
  }
  for(i=L-1; i>0; i--)
  {
    g[i] = g[i-1] + K[i]*f[i-1];
  }
  
  /* filter output */
  return f[L-1];
}


int main(void)
{
  float inp = 1.0f;
  int i;
  FILE* prot = fopen("x.txt", "w");
  
  for(i = 0; i < 100; i++)
  {
    fprintf(prot, "%15f %15f\n", inp, lattice_fir(inp));
    
    inp = 0;
  }
  
  fclose(prot);
  
  return 0;
}

[Johannes E. (cpt_nemo)]

Bert schrieb:
> Vorerst würde ich mich mit einem Lattice FIR begnügen:
> http://de.wikipedia.org/wiki/Latticefilter

Ich sehe jetzt irgendwie nicht, was dein Code mit der Filter-Struktur 
bei Wikipedia gemeinsam hat...

Du berechnest U_old, für was brauchst du das? Und was sind die Speicher 
f und g?

Es gibt bei Wikipedia pro Filterstufe nur einen Speicher (z^-1). Dafür 
braucht man zwei Koeffizienten pro Filterstufe Alpha und Alpha'; bei dir 
gibt es nur ein Koeffizienten-Array K.

[Bert (Gast)]

Hallo Johannes,

ja, ich habe hier im Buch gelesen:

Proakis, J. G. and Manolakis, D. G.
Digital Signal Processing -- Principles, Algorithms and Applications
Pearson International


dass beim Lattice die Koeffizienten alpha und alpha' immer gleich 
sind.
Anscheinend ist eben gerade das der Witz des Lattice; da die Zahlen, die 
man an den Summationspunkten dann addiert, aufgrund von alpha=alpha' 
immer in der selben Grössenordnung sind, ist das ganze numerisch etwas 
besser. (Steht im Buch, lasse ich mal so im Raum stehen, kann ich nicht 
so genau nachvollziehen...)

Ich muss da nochmal drüber schauen, wie man das Filter denn realisieren 
könnte. Oder hast du eine Idee.

Gruss!

[Bert (Gast)]

Hallo,

so nun habe ich auch den Lattice noch geschafft. Ist eigentlich einfach 
:-)
Hier der Code, falls es jemand gebrauchen kann.

Das Lattice lässt sich tatsächlich effizienter implementieren, denn man 
muss nie die alten Messwerte in dem Array herumschieben. Ich denke, das 
ist mit ein Grund, warum Lattice in den Lehrbüchern als so toll 
angepriesen wird.

#define L 4
static const float K[] = {0.1, 0.1, 0.1, 0.1}; /* koeffizienten in matlab berechnen mit tf2latc */

static float g[L];

static void lattice_fir(float input, float* f_out, float* g_out)
{
  float f[L];
  int i;
  
  f[0] = input + K[0]*g[0];
  for(i = 1; i < L; i++)
  {
    f[i] = f[i-1] + K[i]*g[i];
  }
  
  *fout = f[L-1];
  *gout = g[L-1] + K[L-1]*f[L-2];
  
  for(i = L-1; i > 1; i--)
  {
    g[i] = g[i-1] + K[i]*f[i-2];
  }
  
  g[1] = g[0] + K[0]*input;
  g[0] = input;
}