Embedded Files
Fitten: Funktionen an Daten anpassen
Fitten: Funktionen an Daten anpassen
Im letzten Abschnitt haben wir die Daten aus der Datei "scripts/examples/amplitude.dat" in die Tabelle data() geladen (falls das bei Dir nicht so sein sollte, wiederholst Du am besten den letzten Abschnitt). Wir haben außerdem gesehen, dass die Daten Amplituden vs. Frequenz-Kurven enthalten, die jeweils ein Maximum besitzen. Es liegt also nahe, dass es sich hier um Resonanz-Messungen handelt. Als erste Analyse kannst Du also versuchen, eine Resonanzkurve (auch Lorentz-Funktion genannt) an die Daten anzupassen, um die Parameter der Messung zu extrahieren.
Die hier gezeigten Kommandos haben noch viele weitere Optionen, um sie auf spezielle Situationen anzupassen. Das alles hier zu zeigen, würde den Rahmen und die Möglichkeiten sprengen. Wir empfehlen Dir, zu allen Kommandos auch einmal die interne Dokumentation (zum Beispiel help fit und help fft) durchzulesen, um Dich damit vertraut zu machen.
Dazu definieren wir uns zunächst die Lorentz-Funktion. Das ist zwar nicht zwingend nötig, macht den Code im Folgenden aber insgesamt lesbarer:
<- define lorentz(x, x0, amplitude, damping, offset) := amplitude/sqrt((x^2-x0^2)^2+4*damping^2*x^2) + offset
-> Die Funktionsdefinitionen wurden erfolgreich aktualisiert.
Die Dir derzeit noch unbekannten Parameter der Lorentz-Funktion x0, amplitude, damping und offset kannst Du Dir nun mittels des Kommandos fit aus den Daten extrahieren:
<- fit data(:,1:2) -with=lorentz(x, x0, amp, damp, off)
-> Fitte "amp/sqrt((x^2-x0^2)^2+4*damp^2*x^2) + off" ... Erfolg.
===========================================================================================
-> NUMERE: FITERGEBNIS
===========================================================================================
-> Funktion: 0.43003/sqrt((x^2-0.00083679^2)^2+4*3.0264^2*x^2) + 0.0080817
[...]
-> Parameter Initialwert Anpassung Asymptotischer Standardfehler
----------------------------------------------------------------------------------------
amp 0 0.4300348 ± 369.36 (8.589e+004%)
damp 0 3.026374 ± 3.0706e+006 (1.015e+008%)
off 0 0.008081738 ± 4.5169 (5.589e+004%)
x0 0 0.0008367877 ± 2.2077e+010 (2.638e+015%)
[...]
===========================================================================================
Bei Dir wird im Kommandofenster hier noch viel mehr Text stehen. Wir haben die sehr lange Ausgabe, die fit erzeugt, an dieser Stelle etwas abgekürzt und wollen uns erst einmal auf das Ergebnis fokussieren. In der Tabelle findest Du die alphabetisch sortierte Liste der Parameter mit deren Initial- und angepassten Werten sowie einer Schätzung für deren Genauigkeit. Wenn Du die prozentuale Abweichung ansiehst, dann wirst Du sicher vermuten, dass hier etwas nicht stimmt. Falls Du es eher gewohnt bist, mit dem Bestimmtheitsmaß R² zu arbeiten, kannst Du Dir dieses mit der folgenden Zeile auch berechnen lassen:
<- 1 - chi^2 / sum((data(:, 2) - avg(data(:, 2)))^2)
-> ans = 0.08690886
Auch hier zeigt sich, dass offenbar kein besonders guter Match gegeben ist. Die Variable chi wurde dabei automatisch von fit erzeugt. Zur verbesserten Analyse stellen wir einfach mal die Daten mit dem Fit zusammen graphisch dar.
Um das zu tun, verwendest Du wieder das Kommando plot:
<- plot data(), Fit(x)
Es fallen wieder zwei Dinge auf: Du hast data() als Abkürzung für data(:, 1:2) verwenden können (das geht tatsächlich sehr häufig und NumeRe wird diese Abkürzung stets möglichst passend interpretieren). Außerdem hast Du nicht die Funktion lorentz() zum Plotten verwendet, sondern Fit(x). Diese Funktion wird von fit immer mit dem aktuellen Ergebnis erzeugt, so dass man das Ergebnis unkompliziert und schnell darstellen kann.
Trotzdem sieht das Ergebnis aber sehr unschön aus, wie Du auf der rechten Seite erkennen kannst. Zwar wurde eine Funktion durch die Datenpunkte gelegt, aber seltsamerweise scheint kaum ein Punkt wirklich in der Nähe der Funktion zu liegen.
Leider ist das ein fundamentales Problem von Fitting-Algorithmen. Die Ergebnisse sind von den Startwerten abhängig und können durch leichte Änderungen dieser Parameter vollkommen anders ausfallen. Daher ist es sinnvoll, bei nicht-polynomialen Funktionen die ungefähren Parameter als Schätzungen vorzugeben.
Bis jetzt haben wir für alle Parameter beim Fitten den (automatischen) Startwert 0 verwendet (Korrekterweise ist das der automatische Initialwert der neu erzeugten Variablen). Das kannst Du aber auch explizit anders vorgeben, indem Du die params=PARAMS Option verwendest:
<- fit data(:,1:2) -with=lorentz(x, x0, amp, damp, off) params=[x0=2.5,amp=2,damp=0,off=0]
-> Fitte "amp/sqrt((x^2-x0^2)^2+4*damp^2*x^2) + off" ... Erfolg.
===========================================================================================
-> NUMERE: FITERGEBNIS
===========================================================================================
-> Funktion: 0.031963/sqrt((x^2-2.4575^2)^2+4*0.11405^2*x^2) + 0.007953
[...]
-> Parameter Initialwert Anpassung Asymptotischer Standardfehler
----------------------------------------------------------------------------------------
amp 2 0.0319629 ± 0.0087816 (27.47%)
damp 0 0.114046 ± 0.030319 (26.58%)
off 0 0.00795303 ± 0.0045607 (57.35%)
x0 2.5 2.457466 ± 0.015739 (0.6404%)
[...]
===========================================================================================
Wenn Du wieder das Bestimmtheitsmaß R² berechnest, so erhältst Du das folgende:
<- 1 - chi^2 / sum((data(:, 2) - avg(data(:, 2)))^2)
-> ans = 0.9443302
Es werden nun ganz andere Parameterwerte als Ergebnis erzeugt und auch die prozentualen Abweichungen der Fehler sind deutlich kleiner. Wiederholst Du nun den letzten Plot (einfach zweimal [PFEILTASTE-HOCH] und danach [ENTER]), dann wirst Du das Ergebnis auf der rechten Seite bekommen.
Hier liegen die gemessenen Punkte deutlicher in der Nähe der Lorentz-Funktion. Dieser Parametersatz beschreibt die Messpunkte also deutlich besser. Insbesondere die ungefähre Angabe des x0 Parameters hat in diesem Fall eine große Auswirkung auf das Ergebnis. Der Grund hierfür ist, dass x0 im Nenner der Funktion auftritt. Solche Parameter sind besonders empfindlich, da hier kleine Änderungen große Auswirkungen haben können.
Du kannst natürlich auch das Ergebnis für x0 direkt in der Legende anzeigen lassen. Gebe dafür einfach folgendes ein:
<- plot data(), Fit(x) "x_0 = " + #x0
Fouriertransformationen mittels FFT
Fouriertransformationen mittels FFT
Die Frequenz- bzw. Fourieranalyse gehört mit zu dem wichtigsten Handwerkszeug in der Datenanalyse. Daher wollen wir Dir hier auch ein Beispiel zeigen, wie Du eine Fouriertransformation in NumeRe rechnen kannst. Dazu verwenden wir im folgenden das Kommando fft. Zuvor brauchen wir aber noch einen passenden Datensatz, den wir auch transformieren können. Gebe dazu die folgenden vier Zeilen in das Kommandofenster ein:
<- new signal();
<- signal(:, 1) = {0:0.001:0.5};
<- signal(:, 2) = sin(_2pi*50*signal(:, 1)) + sin(_2pi*120*signal(:, 1)) + gauss(0, 2);
<- signal(#, :) = "Time t [s]", "Signal x(t)";
Mit den ersten drei Zeilen haben wir einen etwas verrauschten Datensatz erzeugt, wie er Dir bei der Arbeit mit Audioaufnahmen begegnen kann. In der zweiten Zeile siehst Du außerdem eine Möglichkeit, einen Vektor mit iterativ befüllten Werten zu generieren: {A : B : C} = {A, A+B, A+2*B, ..., C}, wobei C nur dann Teil des Vektors wird, wenn es durch A+n*B ausdrückbar ist. Dabei ist außerdem {A : B} = {A : 1 : B}, bzw. {A : -1 : B}, falls B < A.
Für das Rauschen ist die Funktion gauss(x0,x1) verantwortlich, die Gauss-verteiltes Rauschen um die Position x0 mit der Halbwertsbreite x1 erzeugt. Wir haben außerdem die automatische Vektorisierung durch die Verwendung von Tabellenspalten (signal(:, 1)) genutzt. Die letzte Zeile ergänzt dann, wie Dir bereits bekannt ist, die Spaltenüberschriften.
Stellst Du dann die Daten graphisch dar, zum Beispiel mittels
<- plot signal() -set title="Noisy time domain signal" xlabel="Time t [s]" ylabel="x(t)" interpolate
dann wirst Du zum Beispiel (die Daten sind durch das Gauss-Rauschen zufallsbasiert) das Ergebnis auf der rechten Seite sehen. Dabei haben wir auch zum ersten Mal die Option interpolate verwendet, um die Datenpunkte als eine geschlossene Kurve darzustellen.
Anhand dieser Darstellung auf die in diesem Signal vorhandenen Frequenzen zu schließen, kann recht kompliziert wenn nicht unmöglich sein. Tatsächlich würde man bei diesem Verlauf erwarten, dass es sich nur um Rauschen handelt und jede weitere Analyse unnötig ist.
Dir ist (als erfahrener Data Scientist) aber bekannt, dass man die Frequenzen eines Signals mithilfe einer FFT wieder rekonstruieren kann. Dazu verwendest Du, wie oben bereits angedeutet, das Kommando fft:
<- fft signal()
Das Ergebnis der FFT wird dann in der Tabelle fftdata() angelegt. Du kannst natürlich auch hier die Zieltabelle mit der -target=TABLE() Option angeben. Zur Verbesserung des Signal-To-Noise-Ratios konvertieren wir das Ergebnis dann noch in die sogenannte Power Spectral Density:
<- fftdata(:, 2) ^= 2;
Korrekterweise müsste es fftdata(:, 2) *= conj(fftdata(:, 2)) heißen, aber fft erzeugt standardmäßig eine rein reellwertige Amplituden-Phasen-Repräsentation. Willst Du komplexwertige Ergebnisse haben, musst Du die Option -complex bei fft verwenden.
Stellst Du das Ergebnis dann graphisch dar, kommt es vielleicht zu einer Überraschung:
<- plot fftdata(:, 1:2) -set title="Power spectral density" xlabel="Frequency f [Hz]" ylabel="Power"
Das Spektrum ist horizontal symmetrisch mit positiven und negativen Frequenzen. Außerdem kann man eine horizontale Linie erkennen die von links nach rechts verläuft. Woran liegt das?
Wie Dir vermutlich bekannt ist, enthält ein Frequenzspektrum immer positive und negative Frequenzen, wobei bei einem rein reellen Eingangssignal die negativen Frequenzen exakte Spiegelungen der positiven Frequenzen sind. Die horizontale Linie (bei Dir kann sie auf einer anderen Höhe sein) ergibt sich durch die Tatsache, dass die Frequenzen nicht wie erwartet im Intervall [-500 .. 500], sondern in der Standarddarstellung im Intervall [0 .. 500, -500 .. 0] dargestellt sind und beim Plotten mittels interpolate die 500 und -500-Punkte horizontal (beide Amplituden sind ja identisch) verbunden werden.
Schränken wir also den Zeilenbereich aus der fftdata()-Tabelle ein und verwenden außerdem noch die bars Option, um auf ein Balkendiagramm umzuschalten:
<- plot fftdata(:250, 1:2) -set title="Power spectral density" xlabel="Frequency f [Hz]" ylabel="Power" bars
Wir erkennen nun deutlich zwei scharfe Peaks, die über den Untergrund hinausragen. Beim Anschauen erkennt man, dass deren Frequenzen bei etwa 50 Hz und 120 Hz sind. Um das noch genauer zu verifizieren, können wir die Indices der Tabelle suchen, an denen Amplituden größer als 0.5 sind:
<- logtoidx(fftdata(:250, 2) > 0.5)
-> ans = { 26, 61}
Die Funktion logtoidx() konvertiert dabei die übergebenen Logikwerte (Amplitude größer als 0.5: ja/nein) in die Indices bei denen die logische Aussage wahr ist. Verwenden wir diese beiden Indices, um die zugehörigen Frequenzen zu extrahieren, erhalten wir die folgenden Frequenzen.
<- fftdata({26, 61}, 1)
-> ans = { 50, 120}
Wie wir also vermutet hatten, sind die Peaks bei den Frequenzen 50 Hz und 120 Hz. Da wir die Daten aber genau mit diesen beiden Frequenzen erzeugt hatten, überrascht dieses Ergebnis vermutlich niemanden.
Page updated
Google Sites
Report abuse