Messwerte mit Splines glätten?

Walli

Guten Mittag,

ich hoffe jemand hier hat ein wenig Ahnung davon. Es geht darum recht ungenaue Messwerte durch eine Ausgleichskurve anzunähern. Und zwar hat mein Vater Tabellen mit Messwerten die ungefähr so liegen könnten...

^
|
|
2                                  x
|           x                 x
1   x                  x
|
0---------------1----------------2------>

Er will mit Hilfe dieser Messwerte eine Kurve annähern, um auch z.B. den Wert an der Stelle 1 annähernd bestimmen zu können. Nun habe ich gesagt, dass das kein Problem ist, ich würde das irgendwie hinbekommen. Das war vielleicht ein wenig zu vorschnell!

Mein erster Versuch: Ich habe mit dem Newton-Verfahren das eindeutige Interpolationspolynom bestimmt und Zwischenwerte berechnet. Diese wichen aber zu stark von den erwarteten Werten ab.

Nun habe ich schon einmal etwas von B-Splines gehört. Die sollen doch angeblich keine genau Interpolation sein, sondern eher so etwas wie eine abschnittsweise angenäherte Kurve, die pro Abschnitt mit 4 Punkten konstuiert wird. Ich habe mit Hilfe einer Algorithmensammlung einen Code zusammengebastelt der aber irgendwie nicht wirklich funktioniert. Ich poste mal den relevanten Ausschnitt. Die Stützpunkte stehen in einer deque namens sup, die Elemente der Klasse Point beinhaltet...

// Auswertung an der Stelle x
	double BSpline::eval(const double& x) const {
		// TODO: hier später noch Fehler abfangen

		// suchen, hinter welcher Stützstelle der Wert liegt
		std::deque<Point>::const_iterator it = sup.begin(), end = sup.end();
		while(it->x < x && it != end) ++it;

		// Zur besseren Lesbarkeit Referenzen auf die Punkte holen
		const Point &p0 = *(it - 1), &p1 = *it, &p2 = *(it + 1), &p3 = *(it + 2);

		// Berechnung der Koeffizienten
		double a3 = (-p0.y + 3. * (p1.y - p2.y) + p3.y) / 6.,
			a2 = (p0.y - 2. * p1.y + p2.y) / 2.,
			a1 = (p2.y - p0.y) / 2.,
			a0 = (p0.y + 4. * p1.y + p2.y) / 6.;

		// Wert auf dieses Intervall normieren
		double xrel = (x - p1.x) / (p2.x - p1.x);

		// berechneten Wert herausgeben
		return (((a3 * xrel + a2) * xrel + a1) * xrel + a0);
	}

Btw: Ich weiß selber, dass der Code unter Umständen Bereichsüber-/unterschreitungen und Divisionen durch 0 fabrizieren kann. Das fange ich dann ab wenn es läuft. Im Moment sind die Eingabewerte so konstriert, dass es keine Probleme damit gibt. Ist also ein rein mathematisches Problem hoffe ich ;).

Und nun zeige ich noch ein Bild der Ausgabe...
Ausgabe

Vielleicht hat ja jemand von euch schon einmal etwas ähnliches gemacht und sieht den Fehler. Ich habe leider bisher noch keine Splines in der Numerik-Vorlesung durchgenommen . Alle Informationen über Splines die ich über google etc. finde sind ungefähr so leicht zu verinnerlichen wie ein Telefonbuch. Ich habe auch schon nach fertigen Codes gesucht, aber die sind immer irgendein total unübersichtliches und seitenlanges C.

Gruß,

Christian

Winn

Servus Mastah,

unter -> http://lib-www.lanl.gov/numerical/bookcpdf.html findest Du schonmal SourceCode für Interpolationen. Um die bestmögliche Interpolationsmethode zu benennen, müßte man noch etwas mehr über die Meßwerte wissen.

Werden die Meßwerte stets in gleichen Abständen aufgenommen ? z.B. stets nach der selben Zeitdauer ?
Die Funktionswerte die interpoliert werden sollen, sind jene dann auch immer in den gleichen Abständen ?

Wenn Du beides mit Nein beantwortest, bist du bei den Splines schon sehr richtig. Allerdings sind es dann in der Regel keine B-Splines... sondern einfach nur kubische Splines. Wenn ich mich nicht irre, werden die B-Splines bei der Berechnung von 3D Objekten benutzt. Da schießt Du mit Kanonen auf Spatzen. Die kubischen Splines sind auch in dem obigen Online Book enthalten (Kapitel 3.3 -> Cubic Spline Interpolation )

Gruß Winn

Walli

Hi Winn,

also wie gesagt bin ich noch nicht wirklich mit der Materie vertraut. Ich kann eigentlich nur Interpolieren (Lagrange, Newton, Neville-Aitken), aber das bringt zu schlechte Ergebnisse. Bin mir auch nicht wirklich sicher, ob das B-Splines sind die ich da interpolieren wollte. Da stand irgendwas von Kurvenanpassung und Splines...

Naja, also zu den Messwerten. Die sind in unregelmäßigen Abständen und liegen auf einer relativ glatten Kurve (also ein wenig wellig, aber ohne große Steigungen). Es soll wahllos ein Wert auf dem Messintervall berechnet werden können. Wenn ich eval(x) sage dann will ich auch nur den Wert bei x, und nicht ein ganzes Array von äquidistanten Auswertungen.

Gruß,

Christian

Walli

Also ich habe mir mal das Kapitel zu den kubischen Splines ausgedruckt. Mal schauen, ob das wasfür mich ist...

Danke!

Winn

Servus nochmal,

also dann kommen wirklich nur die Splines für Dich in Frage. Hab das Kapitel nochmal kurz überflogen und muß sagen, auf den paar Seiten versteht man kaum, was die von einem haben wollen. Hier ist mal ein wenig Mathematica Code mit dem ich Splines mal nachgebildet hatte.

Remove["Global`*"]
X = {x0, x1, x2, x3}; // X Werte
Y = {y0, y1, y2, y3}; // Deine Meßwerte
n = Length[X] - 1;
xs[i_] := X[[i + 1]];
ys[i_] := Y[[i + 1]];
s[i_, x_] := a[i]*(x - xs[i])^3 + b[i]*(
x - xs[i])^2 + c[i]*(x - xs[i]) + d[i] // kubischer Spline (da x^3)

g1 = Table[s[i, xs[i]] == ys[i], {i, 0, n - 1}]; // erstes Gleichungssystem, linke Funktionswerte sollen übereinstimmen -> Stetigkeit
g2 = Table[s[i, xs[i + 1]] == ys[i + 1], {i, 0, n - 1}]; // zweites Gleichungsstem, rechte Funktionswerte sollen übereinstimmen -> Stetigkeit
g3 = Table[( D[s[i, x], x] == D[s[i + 1, x], x]) /. 
    x -> xs[i + 1], {i, 0, n - 2}]; // drittes Gleichungssystem, die erste Ableitungen sollen an den Meßwerten übereinstimmen -> Differenzierbarkeit
g4 = Table[( D[s[i,
 x], {x, 2}] == D[s[i + 1, x], {x, 2}]) /. x -> xs[i + 1], {i, 0, n - 2}]; // viertes Gleichungssystem, die zweite Ableitung soll Null sein, dort entstehen dann Wendepunkte
g5 = {(D[s[0, x], {x, 2}] /. x -> 
    xs[0]) == 0, (D[s[n - 1, x], {x, 2}] /. x -> xs[n]) == 0}; // fünftes Gleichungssystem, die Ränder sprich x0 und x3 müssen besonders berücksichtigt werden

L = Solve[{g1, g2, g3, g4, g5} // Flatten, Table[{a[i], 
    b[i], c[i], d[i]}, {i, 0, n - 1}] // Flatten] // Lösen des Gleichungssystems 
sp[x_] := Table[s[i, x] /. L[[1]], {i, 0, n - 1}]

sp[x_] wäre in diesem Fall Dein eval(x), ist aber dort noch in analytischer Form, d.h. erst durch Einsetzen von Meßwerten und ihren entsprechenden x-Werten, entsteht dann Dein kubischer Spline oder auch anders, ein Polynom dritten Grades, welches zufällig (gg) durch die Funktionswerte wandert. Vielleicht wird es so klarer, was da verlangt wird

Gruß Winn

P.S.:
Mathematica Source (mit den Gleichungssytemen) -> http://www.isis.de/~wbilgic/spline.rar
Mathematica Reader -> http://www.wolfram.com/products/mathreader/

Walli

Ups, das sieht ziemlich umfangreich aus. Mal sehen wie ich mir mit den Informationen eine Klasse basteln kann.

Amtrak

hi,
wenn du nur an der Lösung dieses Problems interessiert bist und nicht daran, wie du die Sache mit eigenem Code erschlagen kannst, dann schau dir mal das Programm Gnuplot (unter Linux) an. Was du machen möchtest nennt sich "fitten" von Messwerten. Bei Gnuplot kannst du z.B. ein Polynom x Grades vorgeben, und das Programm versucht dann dieses deinen Messwerten anzunähern. Zum Schluss bekommst du für das Polynom Koeffizienten raus, mit denen du dann Werte an beliebigen Stellen berechnen kannst.

Gruss Alexander

PS: wenns sein muss, habe ich vielleicht sogar noch ein Beispiel, habe ich allerdings nicht direkt zur hand so das ich es suchen müsste. Falls du Interesse hast schau ich mal.

aaaa3

Für Dein Problem sind Splines wahrscheinlich zu kompliziert. In Deinem Fall würde ich mal nach "Approximation" bzw. "Methode der kleinsten Quadrate" googeln.
Du kannst damit ein Polynom bestimmen, dass sich möglichst genau an alle Deine Messwerte annähert. Die Messwerte müssen dabei nicht unbedingt exakt geschnitten werden. Wenn Dei Bild ungefähr die Messwerte widerspiegelt, würde ich vorschlagen, Du versuchst ein kubisches Polynom, oder eine gerade durch Approximation zu bestimmen. Das kubische Polynom wäre etwas genauer, vielleicht reicht Dir aber schon die Gerade.

Walli

Amtrak schrieb:

wenn du nur an der Lösung dieses Problems interessiert bist und nicht daran, wie du die Sache mit eigenem Code erschlagen kannst, dann schau dir mal das Programm Gnuplot (unter Linux) an.

Ich muss da ne Klasse draus machen, da die Werte erst im Programm angegeben werden. Trotzdem danke!

Walli

Noch ne Frage. Muss man bei Splines immer zwingend ein ganzes Array auf einmal berechnen. Geht es nicht, dass man sich 4 Punkte sucht die um den zu interpoliernde Stelle rum liegen und dann mit einer Rechnung auswertet. Weil komischerweise funktioniert meine Lösung mit nem vorberechneten Array und die ohne nicht.

Wenn es die Messwerte sind macht es überhaupt keinen Sinn mit Splines zu interpolieren weil sie
1. rauschen
2. willkürlich sind
3. keiner math. Formel unterliegen

Es gibt aber Möglichkeiten zu interpolieren. Dein Papi muss es dir aber sagen.

Walli

Ich will eben nicht genau interpolieren, das wäre ja kein Problem mit dem Newton-Verfahren. Die Werte bei einer Polynominterpolation weichen zu stark ab, deswegen möchte ich Splines verwenden.

Winn

MaSTaH schrieb:

Noch ne Frage. Muss man bei Splines immer zwingend ein ganzes Array auf einmal berechnen. Geht es nicht, dass man sich 4 Punkte sucht die um den zu interpoliernde Stelle rum liegen und dann mit einer Rechnung auswertet. Weil komischerweise funktioniert meine Lösung mit nem vorberechneten Array und die ohne nicht.

Klar geht das, die Lösung die in dem Mathematica File enthalten ist, ist eine analytische Lösung, d.h. durch Einsetzen Deiner X-Werte und ihren zugehörigen Funktionswerten erhälts Du jeweils abschnittsweise ein kubisches Polynom. Bei vier Stützstellen erhält's Du also drei kubische Polynome, natürlich können es auch nur drei Stützstellen sein, dann haste aber nur zwei Polynome. Ich an Deiner Stelle würde die analytische Lösung innerhalb Deiner Klasse verwenden, einsetzen ist dadurch sehr simpel und der Ablauf dürfte auch einigermaßen flott sein.

Wenn es die Messwerte sind macht es überhaupt keinen Sinn mit Splines zu interpolieren weil sie
1. rauschen
2. willkürlich sind
3. keiner math. Formel unterliegen

Rauschen wird jede Art von Interpolation, die Willkürlichkeit liegt im Rauschen begründet und einer math. Formel unterliegt keiner Interpolation, sondern einer math. Methode ! Interpolationsverfahren, wie Newton, Lagrange steigen rapide im Aufwand je mehr Stützstellen involviert werden. Splines hingegen beruhen auf der Lösung von Gleichungssystemen, die auf nicht vollständig besetzten Matrizen basieren. Dort explodiert der Aufwand nicht gleich.

Walli

aaaa3 schrieb:

Für Dein Problem sind Splines wahrscheinlich zu kompliziert. In Deinem Fall würde ich mal nach "Approximation" bzw. "Methode der kleinsten Quadrate" googeln.

Ups, hab deinen Eintrag gestern überlesen. Sorry .

Ich habe jetzt einen Quellcode zusammengebastelt der ganz gute Resultate bringt. Wenn die Werte nicht zu weit auseinander liegen schneidet die Kurve die Punkte und ansonsten gibt sie eine gute Approximation...

Danke euch allen.

b7f7

man sollte auch versuchen nicht nur potenzen von X im polynom zu verwenden manchmal hat man mit a(1)*x+...+a(n)*cos(x)+a(n+1)*sin(x) bessere erfolge

WebFritzi

Also, ich wollte noch mal sagen, dass Polynome zum Interpolieren meistens scheiße sind. Kubische Splines sind super.

@MastaH: Auf meiner Homepage findest du auch einen kleinen Artikel über kubische Splines.

b7f7

sind kubische spines nicht nur lokale Aproximationen mit polynomen 3.grades?

das Problem mit Polynomen ist eigentlich, dass ein (Rundungs/Rechen)Fehler in den Koeffizienten vor x^25 stärkere wirkung hat als ein Fehler im Linerarglied.
Dies ist ein drung warum man nicht zwischen 2 Punkten mit y=a(1)*x^200+a(0)*x approximiert.

Helium

nicht kubisch, funzt trotzdem:

double interpol (double p0, double p1, double p2, double p3, double p4, double p5, double x)
{
return p2 + 0.04166666666*x*((p3-p1)*16.0+(p0-p4)*2.0
    + x *((p3+p1)*16.0-p0-p2*30.0- p4
    + x *(p3*66.0-p2*70.0-p4*33.0+p1*39.0+ p5*7.0- p0*9.0
    + x *( p2*126.0-p3*124.0+p4*61.0-p1*64.0- p5*12.0+p0*13.0
    + x *((p3-p2)*50.0+(p1-p4)*25.0+(p5-p0)*5.0)))));
}

Winn

WebFritzi schrieb:

Also, ich wollte noch mal sagen, dass Polynome zum Interpolieren meistens scheiße sind. Kubische Splines sind super.

Woow, wir sind einer Meinung ^^

Helium

Für denkfaule:

double hermite (double y0, double y1, double y2, double y3, double x)
{
   return (((1.5f * (y1 - y2) + 0.5f * (y3 - y0)) * x 
      + (y0 - 2.5f * y1 + 2.f * y2 - 0.5f * y3)) * x 
      + 0.5f * (y2 - y0)) * x 
      + y1;
}