Filtern ohne Filter: Vom Oversampling und Tiefpässen

Wer eine Reihe von Meß- oder sonstigen Daten aufgenommen hat und wissen möchte, was zwischen den einzelnen Meßpunkten passiert, muß interpolieren. Eine Möglichkeit, dies zu bewerkstelligen, wird in diesem Beitrag beschrieben. Und wer darüber hinaus schon immer mal wissen wollte, was sich hinter dem Schlagwort Oversampling verbirgt und warum Nachrichtensprecher keine karierten Jacketts tragen, sollte weiterlesen...

Unter Interpolation versteht man bekanntlich die näherungs weise Approximation einer Funktion, von der man nur Stichproben kennt. Diese Stichproben können stündliche Temperaturmessungen, einzelne Punkte einer Bauelementekennlinie oder Amplitudenwerte eines digital abgetasteten Tonsignals sein. Und da Englisch etwas technischer klingt, nennt man die Interpolation im letzteren Falle Oversampling.

Die alte Leier

Mein Gott, fällt denen denn nichts Neues ein, wird sich mancher fragen, das kenne ich doch alles schon. Newton, Lagrange, Spline - alles alte Hüte. In der Tat ist die Interpolation so alt wie die moderne Mathematik. Schon der alte Newton hat ein geniales, für Computer sehr geeignetes Verfahren entwickelt, das der eine oder andere vielleicht noch aus der Schule kennt. Es basiert auf der fundamentalen Erkenntnis, daß durch n Punkte immer genau ein Polynom n-1 ten Grades hindurchgeht. Das heißt, durch zwei Punkte geht immer eine Gerade, durch drei Punkte eine Parabel usw.Bild 1 zeigt in seiner linken Hälfte ein Polynom 7. Grades, das 8 Punkte verbindet. Aber Polynome sind ja nicht das einzige, was die Welt im Innersten zusammenhält. Kubische Splines würden z.B. den Verlauf der 8 Punkte in Bild 1 viel besser approximieren als das Polynom, aber hier ist der Rechenaufwand höher, und letztlich sind sie ja auch Polynome. Ideal wäre ein Verfahren mit niedrigem Rechenaufwand, das nicht so stark schwingt und auch nicht so zackig aussieht wie ein simpler V_PLINE-Aufruf. Der hier vorgestellte Algorithmus verbraucht wenig Rechenzeit und Speicherplatz und, für Programmierer besonders angenehm, er beansprucht wenig Quelltext.

Grau ist alle Theorie

Wollte ein CD-Player mit Oversampling sein Interpolationsproblem mit Polynomen lösen, hätte er bei 44 kHz Abtastfrequenz und 5 Minuten Titellänge ein Polynom 13199999. Grades auszurechnen. Da das auch für den dicksten Signalprozessor zuviel ist, wendet er ein Verfahren an, das ähnlich dem hier beschriebenen funktioniert.

Wie wir aus der Schule wissen, kann man ein periodisches Signal in sogenannte Harmonische, d.h. Grundschwingung und Oberwellen zerlegen (siehe Bild 2). Die Höhe der senkrechten Pfeile in Bild 2 repräsentiert den Betrag des jeweiligen Fourier-Koeffizienten, der seinerseits die Amplitude der n-ten Harmonischen (d.h. Schwingung mit n-facher Grundfrequenz) verkörpert. So ein Pfeil ist in der Signaltheorie eigentlich ein ganz schmales Rechteck, dessen Breite gegen Null und dessen Höhe dafür gegen unendlich geht, so daß es trotz der unendlichen Schmalheit einen Flächeninhalt hat. Da aber das Format der ST-Computer die Darstellung der unendlichen Höhe nicht zuläßt, gibt man mit der Pfeilspitze nur den Betrag des Flächeninhaltes des Rechtecks an, welches in der Signaltheorie Dirne-Stoß genannt wird.

In Bild 3 gehen wir zu einer weniger anschaulichen, aber dafür wissenschaftlich exakteren Darstellung über. Die erste Frage, die einem bei diesem Bild durch den Kopf schießt, ist: Wozu die negativen Frequenzen? Nun, wissenschaftlich unexakt erklärt kann man sich das folgendermaßen vorstellen. Wenn man auf der Zeitachse in positiver Richtung entlangläuft, zählt man soundsoviel Perioden einer Harmonischen pro Sekunde, ergo positive Frequenz. Wenn man hingegen in umgekehrter Richtung läuft, zählt man die gleiche Anzahl von Perioden pro negativer Sekunde, ergo negative Frequenz. Da aber beide Aussagen zutreffen, trägt man kurzerhand positive und negative Frequenz auf. Aber auch die Amplituden sind sowohl beim Vorwärts- als auch beim Rückwärtslaufen die gleichen, deshalb ist der Wert bei f immer genau so groß wie bei -f.

Wenn wir aus dem Frequenzbild die Zeitfunktion wieder zurückgewinnen wollen, zählen wir einfach den Wert von +f und -f zusammen, um auf den Wert der Amplitude der jeweiligen Harmonischen zu kommen. Und da +/und -/vollkommen gleichberechtigt sind, haben sie jeweils die Hälfte der Amplitude der Harmonischen als Betrag (sie sind halb so groß, wie ihr Äquivalent in Bild 2).

Bild 1

Nu gugge an

Die gepunktete Hüllkurve in Bild 3 stellt die Funktion sin(x)/x dar. Sie heißt Spaltfunktion und wird durch si(x) abgekürzt. Wie durch ein Wunder liegen alle Fourier-Darstellungen von Rechteckfolgen unabhängig vom Tastverhältnis in dieser Hüllkurve, es ändert sich nur der Abstand der Dirac-Stöße, wenn man die Breite des Einzelimpulses ti konstant hält und nur die Periodendauer T der Rechteckschwingung ändert. Aber ein Wunder ist das natürlich nicht, sondern Mathematik. Die Hüllkurve ist nämlich das Frequenzspektrum eines einzelnen Impulses aus der Rechteckfolge. Ja richtig, auch nichtperiodische Signale haben ein Frequenzspektrum, wie sollte man sonst einen Rechteckimuls im Telefon hören können, das doch nur Frequenzen von 300 Hz bis 3,4 kHz überträgt. Der Unterschied zur Fourier-Reihe besteht darin, daß das Spektrum eines einzelnen Impulses kontinuierlich ist und somit nicht aus einzelnen Dirac-Stößen besteht. Man könnte sich ein nichtperiodisches Signal z.B. als periodisches mit unendlicher Periodendauer vorstellen. Dadurch wird die Folge der Dirac-Stöße unendlich dicht, und sie werden auch unendlich klein.

Wenn ich aber etwas unendlich Hohes wieder unendlich niedrig mache, kommt etwas Endliches heraus, und das ist gerade das Spektrum eines Einzelimpulses (diesmal ohne.Pfeil, weil endlich).

Heureka!

Und hier haben wir auch schon ein fundamentales Grundgesetz der Signaltheorie entdeckt: Periodifizierung eines Impulses im Zeitbereich wird zur Abtastung (und nichts anderes sind Dirac-Stöße, die in einer Hüllkurve liegen) des Spektrums dieses Impulses im Frequenzbereich, wobei der Abstand der Dirac-Stöße genau der Grundfrequenz (1/T) der Periodifizierung entspricht. Wobei wir unter Impuls eine beliebige zeitabhängige Funktion verstehen.

Aber was dem Frequenzbereich recht ist, kann dem Zeitbereich nur billig sein. Denn aus der Mathematik folgt, daß auch Periodifizierung im Frequenzbereich zur Abtastung im Zeitbereich wird, und gerade das ist es ja, was wir in unserem Zeitbereich haben, wenn wir interpolieren: Abtastwerte. Demzufolge müßte das Frequenzspektrum der Abtastwerte eine periodische Folge des Spektrums der Originalfunktion sein (siehe Bild 4), die wir ja mit unserer Interpolation finden wollen. Und das müßte nicht nur so sein, es ist auch so.

Wenn es uns nun also gelänge, das Originalspektrum aus dem periodischen Spektrum herauszufiltern, hätten wir eine eindeutige Repräsentation der Originalfunktion im Frequenzbereich. Doch hier taucht ein Problem auf, das jeder Programmierer von Spielen mit Digisound kennt. Die „Periodendauer“ im Frequenzspektrum ist -wer hätte es anders erwartet - der Kehrwert des Abstandes zweier Abtastwerte im Zeitbereich, so wie die Periodendauer im Zeitbereich ja auch der Kehrwert des Abstandes zweier Dirac-Stöße im Frequenzbereich ist (siehe Bild 3). Somit werden die Abstände der Teilfunktionen in Bild 4 mit Größerwerden der Abtastzeit immer geringer.

Sollte die Abtastzeit zwischen zwei Proben (,Samples') zu groß werden, überlappen sich die einzelnen Teilspektren, und eine eindeutige Identifikation des Spektrums der Originalfunktion ist nicht mehr möglich (siehe Bild 5). Dieses Hineinreichen eines weiter oberhalb liegenden Spektrums in das Originalspektrum hören wir bei schlechtem Digisound als hohes Zwitschern im Hintergrund oder wir sehen es als Moire auf dem karierten Jackett des Nachrichtensprechers, wenn die „Frequenz“ des Musters zu hoch ist für die „Abtastrate“ der Fernsehkamera. Nur daß sich in diesem Falle Frequenz und Abtastrate auf den Ort und nicht auf die Zeit beziehen.

Bild 2

Schreiten wir zur Extraktion...

Nehmen wir einmal an, die Abtastzeit sei kurz genug, so daß sich die Spektren nicht überlappen. Wenn das der Fall ist, kann man das Originalspektrum sehr leicht her-ausfiltern, indem man den ganzen Salat mit einer nichtperiodischen Rechteckfunktion multipliziert (Bild 6). Herausfiltern ist dabei wörtlich gemeint, denn so ein Rechteck im Frequenzbereich ist ja nichts anderes als ein idealer Tiefpaß, der unterhalb einer gewissen Frequenz alles durchläßt und oberhalb dieser Frequenz nichts. Jetzt haben wir das Spektrum der Originalfunktion, wie es leibt und lebt, und brauchen, da es eine eindeutige Abbildung derselben ist, nur noch die Rücktransformation vorzunehmen. Aber gerade hier ist der Haken: Die Rücktransformation ist nämlich sehr aufwendig und würde den Rahmen dieses Artikels sprengen. Doch Murphy wird diesmal überlistet, es gibt nämlich eine Regel, die besagt, daß die Zeitfunktion eines Produktes im Frequenzbereich (zwischen Rechtecktiefpaß und periodischer Spektrenfolge) gleich der Faltung der beiden Zeitfunktionen der Faktoren ist (was Faltung ist, kommt gleich). Das heißt, ich nehme die Zeitfunktion des periodischen Spektrums (unsere Abtastfolge, aus der immer noch nichts gescheites geworden ist) und die Zeitfunktion des Filterrechtecks und falte sie miteinander.

Rüber und nüber

Faltung ist leider nicht das einfache Herumklappen irgendeines Funktionsteils auf einen anderen, sondern auch wieder eine komplizierte Operation, die zu erklären den Rahmen des Artikels sprengen würde. Wir haben aber Glück, da die Faltung einer Funktion (der Zeitfunktion des Tiefpasses) mit einem Dirac-Stoß (einem unserer Meßwerte) nichts weiter ist, als die Verschiebung dieser Funktion an die Stelle des Dirac-Stoßes und Multiplikation derselben mit dem Betrag des Dirac-Stoßes (Bild 7 links). Und das fällt unserem Atari gar nicht schwer. Wenn wir also für alle Meßwerte, respektive Dirac-Stöße, diese Operation ausführen und alles zusammenaddieren (Bild 7 rechts), haben wir die Spektralfunktionen multipliziert, ohne sie jemals gesehen zu haben, und die Originalfunktion liegt vor. Das einzige Problem ist die Kenntnis der Zeitfunktion des Tiefpaßfilters und die ist - Bild 7 hat es schon verraten - auch wieder eine Spaltfunktion (wir erinnern uns an die Spektralfunktion des Rechteckimpulses). Unter der Maßgabe, daß die Abtastzeit zwischen zwei Werten kurz genug und immer einheitlich lang war und sich die periodischen Spektren auch nicht überlappen, haben wir nun wieder die Originalfunktion wie sie vor der Abtastung war (Bild 1 rechts), und sind glücklich und zufrieden.

Wenn der CD-Player also zweifach oversamplet, interpoliert er jeweils einen Wert zwischen zweien, die wirklich auf der CD stehen. Oversampling ist demnach nichts weiter als Tiefpaßfilterung, nur daß die Grenzfrequenz des Tiefpasses über dem Hörbereich liegt. Aber wir sagen lieber weiterhin Oversampling, da Tiefpaßfilterung so werbewirksam ist, wie die Feststellung, daß Edelkäse aus verfaulter Milch gemacht wird.

In medias res

Jetzt wollen wir aber endlich auch mal was programmieren. Die Spulen, Kondensatoren und Widerstände lassen wir im Bastlerladen und konstruieren unsere Filter mit RIEM ANN2.PRG. Wie vielleicht nicht jeder weiß, ist Riemann nicht nur ein hervorragendes Mathematikprogramm, sondern auch eine richtige Programmiersprache, in der der Liste. Somit steht nun in DX der Abstand zwischen zwei Abtastwerten und in B offensichtlich die Anzahl der Meßpunkte. Initial-7 ist Null und in der folgenden Zählschleife in A von 1 bis B wird für jeden Meßpunkt die oben besprochene Verschiebe- und Multiplikationsoperation durchgeführt und aufsummiert. FIRST(NTHNODE( YLIST,A)) bezeichnet das A-te Element von YLIST, also den Betrag des Dirac-Stoßes, und FIRST (NTHNODE(XLIST,A)) das A-te Element von XLIST, also die Stelle, an die wir die Gewichtsfunktion (so heißt nämlich die Zeitfunktion des Tiefpasses) verschieben müssen.

Ergebnis einer jeden Funktion ist immer der zuletzt ausgewertete Ausdruck, also letztlich Y in INTERP. Der ideale Tiefpaß in der Gewichtsfunktion G ist in Kommentarklammern eingeschlossen, da ein Kosinusquadrattiefpaß meist bessere Ergebnisse liefert. Nun aber noch ein Wort zur Spaltfunktion SI(X): Natürlich würde Riemann aufgrund der bereits vorhandenen Intelligenz den Grenzwert sin(0)/0 auch alleine herausbekommen, aber mit der angegebenen Konstruktion geht es schneller. Warum diese Konstruktion identisch zu ,if X=0 then 1 else SIN(X)/X‘ ist, mag der geneigte Leser selbst herausbekommen. Ein kleiner Tip: für die booleschen Operatoren gibt es kein TRUE, sondern nur FALSE oder nicht-FALSE.

Bild 3

Bild 4

Fazit

Auch wenn in diesem Beitrag immer von Zeitbereich und Zeitfunktion die Rede war, kann man statt der Zeit natürlich jede andere Maßeinheit genausogut einsetzen - siehe Fernsehkamera. Digitalfilter in diversen Soundprogrammen, mit denen man die Klangfarbe einstellen kann, arbeiten übrigens nicht nach diesem Verfahren, sie verändern nur die schon vorhandenen Werte mit sogenannten Kammfiltern, welche allerdings ähnlich funktionieren. Wem das alles nicht wissenschaftlich genug war, oder wer einmal andere Tiefpässe oder Bandpässe ausprobieren möchte, dem ist mit [1] sicher weitergeholfen. Und wer [2] besitzt, ist vielleicht gar nicht bis zu dieser Textstelle vorgestoßen. Und nun viel Spaß beim Experimentieren.

Literatur:

[ 1 ] Kreß/Irmer Angewandte Systemtheorie, Verlag Technik 1989
[ 2 ] Riemann II-Handbuch

% Diese Funktion liefert den interpolierten Wert an der 
  Stelle X aus den durch XLIST und YLIST spezifizierten Punkten. %

FUNCTION INTERP(X, XLIST, YLIST,%Lokal:%DX,Y,A,B),
    DX:SECOND(XLIST) - FIRST(XLIST),
    B: LENGTH(XLIST),
    Y:0,
    NLOOP A,1,B,
        Y:Y+FIRST(NTHNODE(YLIST,A))
        *G(X-FIRST(NTHNODE(XLIST,A))),
    ENDNLOOP,
    Y ENDFUN$

% Das ist die Zeitfunktion des Tiefpasses. Durch Aufhebung 
  der Kommentarklammern kann man zwischen einem idealen 
  oder einem weicheren cos*2-Tiefpass wählen. %

FUNCTION G(X),
    % SI(##PI*X/DX)                     IDEALER TIEFPASS %
    SI(2*##PI*X/DX)/(1-(2*X/DX)A2) %    COS^2 - TIEFPASS %
ENDFUN$

% Die Spaltfunktion. %

FUNCTION SI(X),
    X=0 AND 1 OR SIN(X)/X
ENDFUN$

Georg Michel
Links

Copyright-Bestimmungen: siehe Über diese Seite