Fehlerminimierung der Spline-Berechnung auf FPGA

[Christoph K. (odsource)]

Schönen guten Tag alle zusammen,

ich arbeite momentan an der Synchronisation von Datenströmen auf einem 
FPGA (noch simulativ mit REAL) und nutze hierfür die Interpolation 
mittels natürlicher kubischer Splines.
Hierfür erhalte ich von Data Clock Recovery (DCR) Modulen jeweils die 
Daten jedes Datenstroms. Erhalte ich das enable eines DCR speichere ich 
den erhaltenen Wert mit einem Zeitstempel (intern generiert) ab. Sobald 
ich initial 4 Messwerte pro Kanal erhalten habe führe ich die 
Interpolation jeden Takt aus. Ich berechne immer nur das Polynom 
zwischen den mittleren beiden Messwerten.
Dies funktioniert auch soweit, jedoch ist der Fehler vergleichsweise 
groß. Ich habe den VHDL Code in Matlab übertragen und erhalte eine 
Abweichung von ca. +/- 500. Wenn ich jedoch eine nicht lineare 
Regression mit least square (ebenfalls kubisch) in Matlab nehme lande 
ich bei einem Fehler von +/- 4. Auf den selben Fehler komme ich auch mit 
der spline() Funktion in Matlab. Diese nutzt jedoch mehr Koeffizienten 
(obwohl in der Dokumentation steht das es kubisch ist). Hat jemand eine 
Idee wie ich vorgehen kann? Sollte vom Ressourcenverbrauch nicht zu 
überschwänglich sein da ich 100 Takte Zeit habe. Eventuell kann ich es 
auf 200 Takte erweitern aber das wäre eher die Notlösung.

Im Anhang hab ich mal alle Matlab Dateien angehängt.

Mit freundlichen Grüßen
Christoph

[Duke Scarring (Gast)]

Was sind denn das für Datenströme?
Wie schnell kommen die?
Welche Bitbreiten habe sie?
Und was wird als eigentliches Resultat erwartet?

[A. B. (Gast)]

Christoph K. schrieb:
> Dies funktioniert auch soweit, jedoch ist der Fehler vergleichsweise
> groß. Ich habe den VHDL Code in Matlab übertragen und erhalte eine
> Abweichung von ca. +/- 500. Wenn ich jedoch eine nicht lineare
> Regression mit least square (ebenfalls kubisch) in Matlab nehme lande
> ich bei einem Fehler von +/- 4. Auf den selben Fehler komme ich auch mit
> der spline() Funktion in Matlab. Diese nutzt jedoch mehr Koeffizienten

Kubische Spines sind sehr gutmütig, das zugehörige LGS auch. Wenn da 
also unerklärlich große Fehler herauskommen, ist der Code nicht korrekt.
Im Zweifel rechnet man die Koeffizienten mal per Hand nach und 
vergleicht:

0 & 50 & 0.000155933(x-50)^{3} & +0(x-50)^{2} & +12.3307(x-50) & -1845 
\\
1 & 150 & -0.000435667(x-150)^{3} & +0.04678(x-150)^{2} & 
+17.0087(x-150) & -456 \\
2 & 250 & 0.000279733(x-250)^{3} & -0.08392(x-250)^{2} & +13.2947(x-250) 
& +1277 \\

Außerdem sind natürliche kubische Splines eindeutig bestimmt, ebenso der 
Begriff als solcher. Wenn also in der Matlab-Doku ausdrücklich davon die 
Rede ist, kann es keinen Zweifel geben, dass dort auch (im Rahmen der 
Rechengenauigkeit) genau die oben angegebenen Koeffizienten herauskommen 
müssen.

[Christoph K. (odsource)]

Vielen Dank für die schnellen Antworten.

Duke Scarring schrieb:
> Was sind denn das für Datenströme?
> Wie schnell kommen die?
> Welche Bitbreiten habe sie?
> Und was wird als eigentliches Resultat erwartet?

Es soll später quasi wie ein Oszilloskop funktionieren. Also sollen 
prinzipiell alle Arten von Datenströmen möglich sein. In der Simulation 
nutze ich allerdings einen generierten Sinus mit 20kHz. Die Frequenz 
soll allerdings auch variabel sein (im Rahmen des Nyquist Theorems). Ich 
arbeite momentan mit 100 MHz. Es sollen 20 Bit Fix Point Daten sein. 16 
Vor- und 4 Nachkommabits. Als Resultat erwarte ich eine möglichst 
geringe Abweichung zum "echten" Wert. Wenn ich (auf 2 Kanäle gesehen) 
bspw. von Kanal 1 zwischen zwei Messwerten mit den Zeitstempeln 100 und 
200 und von Kanal 2 zwischen den Messwerten mit den Zeitstempeln 130 und 
230 jeweils interpoliere sollte ich ja - so zumindest in der Theorie - 
mit einem x Wert zwischen 130 und 200 einen Zeitpunkt erwischen an dem 
beide Datenströme synchron sind.

A. B. schrieb:
> Kubische Spines sind sehr gutmütig, das zugehörige LGS auch. Wenn da
> also unerklärlich große Fehler herauskommen, ist der Code nicht korrekt.
> Im Zweifel rechnet man die Koeffizienten mal per Hand nach und
> vergleicht

Vielen Dank für den Tipp, ich hatte tatsächlich einen Fehler bei der 
Berechnung von a, jedoch hilft mir das nur in diesem konkreten Beispiel 
auf einen Fehler von ca 8 zu kommen. Wenn ich vom selben Signal die 
x-Werte um 50 Zeiteinheiten schiebe und die jeweiligen y-Werte nehme bin 
ich bei einem Fehler von ca. 800 wenn ich auch für den Zeitpunkt 175 das 
ganze berechne.

x=[50; 150; 250; 350];
y=[-1845; -456; 1277; 2047];

175 liegt ganz genau auf dem Wert 0, deshalb verwende ich immer diesen 
Punkt für die Berechnung.

A. B. schrieb:
> Wenn also in der Matlab-Doku ausdrücklich davon die
> Rede ist, kann es keinen Zweifel geben, dass dort auch (im Rahmen der
> Rechengenauigkeit) genau die oben angegebenen Koeffizienten herauskommen
> müssen.

Wenn ich mir in Matlab via

pp = spline(x,y);
fprintf('coef=%.6e\n',pp.coefs);

 die Koeffizienten ausgeben lasse erhalte ich allerdings all die 
folgenden Koeffizienten:

>> SplineMethod
coef=-2.178333e-04
coef=-2.178333e-04
coef=-2.178333e-04
coef=8.255000e-02
coef=1.720000e-02
coef=-4.815000e-02
coef=7.813333e+00
coef=1.778833e+01
coef=1.469333e+01
coef=-1.845000e+03
coef=-4.560000e+02
coef=1.277000e+03

Nachtrag:
Hab gerade meinen Fehler bei der Denkweise gefunden. Es liefert für alle 
drei Polynome die Koeffizienten zurück :-P
Ich war so auf meine Berechnung zwischen zwei Messwerten fixiert das ich 
es jetzt erst gesehen habe.

[A. B. (Gast)]

Christoph K. schrieb:
> Vielen Dank für den Tipp, ich hatte tatsächlich einen Fehler bei der
> Berechnung von a, jedoch hilft mir das nur in diesem konkreten Beispiel
> auf einen Fehler von ca 8 zu kommen. Wenn ich vom selben Signal die
> x-Werte um 50 Zeiteinheiten schiebe und die jeweiligen y-Werte nehme bin
> ich bei einem Fehler von ca. 800 wenn ich auch für den Zeitpunkt 175 das
> ganze berechne.

Das verstehe ich nicht: Wenn man alle x-Werte verschiebt, ändert sich an 
den Koeffizienten des Splines gar nichts. Die Polynome des Splines sind 
gerade so gebaut, dass immer nur die Differenz von der Variablen x zum 
x-Wert des betreffenden Stützpunkt darin vorkommen. Und auch bei der 
Berechnung sieht man doch, dass nur die Differenzen der x-Werte 
vorkommen.

> coef=-2.178333e-04
> coef=-2.178333e-04
> coef=-2.178333e-04
> coef=8.255000e-02
> coef=1.720000e-02
> coef=-4.815000e-02
> coef=7.813333e+00
> coef=1.778833e+01
> coef=1.469333e+01
> coef=-1.845000e+03
> coef=-4.560000e+02
> coef=1.277000e+03

Die letzten drei stimmen, das sind die konstanten Terme. Aber der Rest 
passt überhaupt nicht. Dann macht Matlab da irgendetwas anderes ... 
Irgendwo müssten doch die natürlichen Randbedingungen angegeben werden?

[Christoph K. (odsource)]

A. B. schrieb:
> Das verstehe ich nicht: Wenn man alle x-Werte verschiebt, ändert sich an
> den Koeffizienten des Splines gar nichts. Die Polynome des Splines sind
> gerade so gebaut, dass immer nur die Differenz von der Variablen x zum
> x-Wert des betreffenden Stützpunkt darin vorkommen. Und auch bei der
> Berechnung sieht man doch, dass nur die Differenzen der x-Werte
> vorkommen.

Sorry ich weiß nicht ob ich es richtig verstehe. Wenn ich das Polynom 
zwischen den x-Werten 150 und 250 berechne bekomme ich doch andere 
Koeffizienten wie wenn ich es zwischen den x-Werten 100 und 200 
berechne, oder?

A. B. schrieb:
> Irgendwo müssten doch die natürlichen
> Randbedingungen angegeben werden?

Matlab nutzt die not-a-knot Bedingung für die Berechnung, jedoch sollte 
das für das "mittlere" Polynom keinen Unterschied machen, sondern nur 
für die anderen beiden. Mich interessiert ja nur das "innere" Polynom.

[A. B. (Gast)]

Christoph K. schrieb:
> Sorry ich weiß nicht ob ich es richtig verstehe. Wenn ich das Polynom
> zwischen den x-Werten 150 und 250 berechne bekomme ich doch andere
> Koeffizienten wie wenn ich es zwischen den x-Werten 100 und 200
> berechne, oder?

Es gibt nicht "das" Polynom, sondern es sind drei Stück. Nein, die 
Koeffizienten hängen ausschließlich von den Differenzen der x-Werte 
untereinander ab. Etwas anderes ist hinterher das Auswerten des Splines 
für bestimmte x-Werte, sprich Einsetzen eines x in eines der vorher 
berechneten Polynome. Klar, wenn man da andere Werte einsetzt oder für 
das betreffende x dann ein anderes der drei "zuständig" wird ...

Eine Verschiebung der x-Werte am Anfang ist ja genauso, als ob man die 
Grafik (s. o.) in horizontaler Richtung verschiebt (bzw. nur die 
Beschriftung an der x-Achse): Am Bild ändert sich nichts. Was man 
hinterher einsetzt, sagt dagegen, wo man in der Zeichung hinschaut.

> Matlab nutzt die not-a-knot Bedingung für die Berechnung, jedoch sollte
> das für das "mittlere" Polynom keinen Unterschied machen, sondern nur
> für die anderen beiden. Mich interessiert ja nur das "innere" Polynom.

Oh doch, dadurch sieht das Gleichungssystem völlig anders aus und es 
kommen dann natürlich völlig andere Polynome heraus. Zwar wird der 
Unterschied weit "außen" natürlich größer sein, aber prinzipiell ist's 
dann überall anders.
Der Witz bei der Spline-Interpolation ist doch gerade, dass die Polynome 
aneinander gekoppelt sind. Ändert man irgendwo etwas, hat das überall 
Auswirkungen. Das Bild mit der Holzlatte: Wenn man die z. B. am Rand 
biegt, verformt sich die überall.

[Christoph K. (odsource)]

A. B. schrieb:
> Es gibt nicht "das" Polynom, sondern es sind drei Stück. Nein, die
> Koeffizienten hängen ausschließlich von den Differenzen der x-Werte
> untereinander ab.

Das drei Polynome rauskommen stimmt, jedoch interessiere ich mich immer 
nur für das Zweite. Im Anhang hab ich mal "händisch" 
("Bsp_Calc_2_Channel.txt") die zwei Beispiele gerechnet und bekomme da 
auch vollkommen andere Koeffizienten (wieder nur für das "mittlere" bzw. 
zweite Polynom gerechnet). Genutzt habe ich die Methode von dem 
angehängten Paper ("spline_note.pdf") mit natürlichen Splines.

> Oh doch, dadurch sieht das Gleichungssystem völlig anders aus und es
> kommen dann natürlich völlig andere Polynome heraus. Zwar wird der
> Unterschied weit "außen" natürlich größer sein, aber prinzipiell ist's
> dann überall anders.

Das stimmt wohl, jedoch bekomme ich auch andere Koeffizienten als Matlab 
wenn ich es von Hand berechne (siehe "Not-a-knot_Spline.txt"). Und das 
Ergebnis ist dem mit den natürlichen Bedingungen sehr ähnlich. Auf dem 
Bild "Spline_Vergleich" hab ich alle drei Varianten mal 
gegenübergestellt. Die linke Spalte ist für die x-Werte 0, 100, 200, 300 
und die linke für 50, 150, 250, 350. Die unterste Spalte jeder 
Berechnung ist das Ergebnis. Rundungsfehler seien mal dahingestellt aber 
das sollte gerade auf der linken Seite nicht so einen großen Unterschied 
machen.

[A. B. (Gast)]

Christoph K. schrieb:
> Das stimmt wohl, jedoch bekomme ich auch andere Koeffizienten als Matlab
> wenn ich es von Hand berechne (siehe "Not-a-knot_Spline.txt"). Und das
> Ergebnis ist dem mit den natürlichen Bedingungen sehr ähnlich. Auf dem
> Bild "Spline_Vergleich" hab ich alle drei Varianten mal
> gegenübergestellt. Die linke Spalte ist für die x-Werte 0, 100, 200, 300
> und die linke für 50, 150, 250, 350. Die unterste Spalte jeder

Hm, ich sehe da vier Spalten mit Koeffizienten, was ist was? Und wo 
kommt da die -1277 her??? Der konstante Koeffizient beim 2. Polynom (für 
150 bis 250) muss immer -456 sein, denn der ist doch gerade der y-Wert 
vom "linken" x-Wert.

[Christoph K. (odsource)]

A. B. schrieb:
> Hm, ich sehe da vier Spalten mit Koeffizienten, was ist was? Und wo
> kommt da die -1277 her??? Der konstante Koeffizient beim 2. Polynom (für
> 150 bis 250) muss immer -456 sein, denn der ist doch gerade der y-Wert
> vom "linken" x-Wert.

Sorry, ganz vergessen das Bild richtig zu erklären. Auf der linken Seite 
stehen die Koeffizienten, in der Mitte der berechnete Wert mit diesem 
Koeffizienten und rechts steht die Potenz. Somit steht links unten der 
Koeffizient a und nach oben hin geht es bis d. Berechnet wird das ganze 
immer für den x-Wert 175. Sowohl in der linken als auch der rechten 
Spalte.

Auf der linken Seite stehen die Koeffizienten für folgende x-y 
Wertepaare

x=[0; 100; 200; 300];
y=[-2047; -1277; 456; 1845];

auf der rechten Seite stehen die Koeffizienten für diese Paare

x=[50; 150; 250; 350];
y=[-1845; -456; 1277; 2047];

Und laut dem Paper berechnet sich der Koeffizient a wie folgt:

$$a_i = y_i$$

Wenn ich also die Polynome für die oberen x-y Wertepaare berechne 
bekomme ich für a beim zweiten Polynom den Wert -1277 und bei den 
anderen Wertepaaren ein a mit dem Wert -456 für das zweite Polynom.

Vielleicht auch nochmal fürs Verständnis. Die linke Spalte stellt einen 
Datenstrom dar und die rechte einen Zweiten. Sie sind somit quasi 
unabhängig voneinander, stellen allerdings das gleiche Signal dar. Somit 
muss ich einmal für den einen und einmal für den anderen Datenstrom das 
zweite Polynom berechnen.

PS: Vielen Dank nochmal bin echt positiv überrascht wie gut ein Forum 
funktionieren kann :-)

[Christoph K. (odsource)]

Für alle die es interessiert. Ich hab meine beiden Fehler gefunden.

1. In beiden Fällen habe ich im Polynom

$$a + b*(x-x_i) + c*(x-x_i)^2 + d*(x-x_i)^3$$

 statt

$$(175 - x_i)$$

 immer mit

$$(175 - 150)$$

, also dem

$$x_i$$

 von nur einer Spline Berechnung gerechnet.

2. Meine Gauss Implementierung ist für die not-a-knot Bedingung zu 
ungenau. Durch Rundungsfehler erhalte ich insgesamt eine Abweichung von 
11 statt 3,7.

Wenigstens weiß ich jetzt das die not-a-knot Bedingung um einiges 
präziser ist als ein natürliches Spline.