Datenanpassung an trigonometrische FunktionenHallo,
Ich kämpfe damit, gemessene Daten an eine trigonometrische Ksin(ax+b)+c-Funktion anzupassen. Dies ist die Natur meines Prozesses und ich suche nach einer automatischen Möglichkeit, ihn durch JSL-Code anzupassen.
Ich habe einen Code geschrieben, der die Daten normalisiert, sodass die Oszillation zwischen -1 und 1 liegt.
Am Ende möchte ich a und b Konstanten erhalten (K und c auch, aber es ist viel einfacher)
angehängtes Beispiel für meinen Datentyp.
BR,
Yoav
AnalytikDatentabelleSonstigesSkriptingRe: Datenanpassung an trigonometrische Funktionen FunctionIch weiß nicht, ob Sie das wollen, aber Sie können dies mit der nichtlinearen Plattform lösen.Angehängt sind Ihre Daten, zusammen mit einer Prädiktorspalte (die die Parameter zum Schätzen verwendet) und einem Skript, das nach diesen Parametern auflöst.Ich gebe zu, dass es eine ganze Weile gedauert hat, erste Vermutungen über die Parameter zu sammeln, die die Konvergenz der Lösung ermöglicht haben (wahrscheinlich aufgrund meiner mangelnden Vertrautheit mit der Sinusfunktion - es ist lange her, dass ich mit solchen Funktionen gearbeitet habe).Aber konzeptionell wird die nichtlineare Plattform das erreichen, was Sie wollen.
Re: Datenanpassung an trigonometrische Funktionen FunctionSo erhalten Sie eine Schätzung der ein Parameter mit der FFT-Funktion. Die Kommentare zu Sekunden und Stunden gehen davon aus, dass die Daten etwa alle 4 Sekunden für etwa eine Stunde erfasst wurden; die FFT kümmert sich nicht um die Einheiten. Die FFT ordnet die Frequenzen in den Daten so, dass das erste Bin nach dem DC-Offset-Bin die niedrigste Frequenz darstellt, die einen vollständigen Zyklus in den Eingabedaten ausmachen kann, das nächste Bin 2 vollständige Zyklen umfasst, ... das mittlere Bin ist ndata/ 2 Zyklen. Die Behälter danach sind ein Spiegel des ersten.
dt = Open( "z:\example data.jmp" );
// by hand, I made a formula column to look at x-lag(x) and
// decided the first 79 and last 8 were not representative...
rowstart = 80;
rowstop = 937;
timeStart = dt:x[rowstart];
timeStop = dt:x[rowstop];
ydata = dt[rowstart :: rowstop, {y}]; // throw away a few with irregular time
{r, i} = FFT( {ydata} );
nBins = N Items( r );
sampleRate = (timeStop - timeStart) / (nBins - 1); // imagine there were only 2 bins...thus -1.
// the samplerate is ~4 seconds/sample...
// the highest frequency this can resolve is 1 cycle/ 8 seconds. (.125 Hz)
// frequencies higher than that will cause aliasing unless they are filtered out before this.
// the bin [nBins/2] represents the maximum frequency that can be resolved
lastBinFrequency = 1 / (2 * sampleRate); // ~ 0.125 Hz
// the first bin [1] is the DC, or offset from 0, component. We'll ignore it.
// the next bin [2] represents the frequency that has one cycle over the length of the data.
// the lowest non-DC frequency is one cycle per hour (the length of time the data was collected).
binWidth = 1 / (timeStop - timeStart); // ~0.0003Hz, or ~1 cycle per hour
// note: lastBinFrequency == (nbins-1)*binWidth/2 -- each bin is this wide, in frequency
// the remaining bins are a mirror of the first and ignored.
// combine real and imaginary bin components to get an amplitude
p = Sqrt( r ^ 2 + i ^ 2 );
dtGraph2 = As Table( ((rowstart :: rowstop) - rowstart)` || p );
dtGraph2 << Graph Builder(
Size( 1300, 500 ),
Show Control Panel( 0 ),
Show Legend( 0 ),
Variables( X( :Col1 ), Y( :Col2 ) ),
Elements( Points( X, Y, Legend( 5 ) ) ),
SendToReport(
Dispatch( {}, "Col1", ScaleBox, {Format( "Best", 12 ), Min( 0 ), Max( nBins / 2 ), Inc( 100 ), Minor Ticks( 4 )} ),
Dispatch( {}, "Col2", ScaleBox, {Format( "Best", 12 ), Min( 0 ), Max( 30000 ), Inc( 5000 ), Minor Ticks( 1 )} ),
Dispatch( {}, "graph title", TextEditBox, {Set Text( "FFT Bins" )} ),
Dispatch( {}, "X title", TextEditBox, {Set Text( "bin number" )} ),
Dispatch( {}, "Y title", TextEditBox, {Set Text( "size" )} )
)
);
// there is a stong DC component in this data (all values > 0) so bin 1 is the max.
If( Loc Max( p ) != 1,
Throw( "something odd" )
);
p[1] = 0; // clear the DC value
// next max after DC:
maxBinIndex = Loc Max( p ); // 5. the data repeats about (5-1)==4 times in this frame.
Write(
"\!nBiggest frequency component at ",
binWidth * (maxBinIndex - 1),
" Hz or ",
Round( 1 / (binWidth * (maxBinIndex - 1)) ),
" seconds per cycle"
);
// where does a==.007 come from?
secondsPerCycle = 1 / (binWidth * (maxBinIndex - 1));
totalSeconds = timeStop-timeStart;
cyclesInSample = totalSeconds/secondsPerCycle;//about 4
timerange = timeStop-timeStart;
a_from_fft = cyclesInSample*2*pi()/timerange; // .0073...
show(a_from_fft);
Eine starke Reaktion im 4. Bin nach DC-Bin
Größte Frequenzkomponente bei 0,00116448326055313 Hz oder 859 Sekunden pro Zyklus
a_from_fft = 0,007316664113164;
Die DC-Komponente bezieht sich auf den Mittelwert der Daten, die C Wert.
Re: Datenanpassung an trigonometrische Funktionen FunctionIch denke, es hilft, das zu verwenden, was ich als Standardform einer Sinusfunktion bezeichne: Y = K*sin(A*(x -
) + C.
In Ihren Daten:
K = 1/2 des Wertebereichs
A = 2*pi / Periode
B = Verschiebung nach rechts von 0 zum Startpunkt der Sinuswelle
C = Mittelwert Ihrer Y-Werte
Ausgangsschätzungen für die nichtlineare Regression können mithilfe dieser Definitionen von der Distributionsplattform und dem Graph Builder gefunden werden:
K = 1/2 * (Max - Min) = 1/2 * (47 - 17,6) = 14,7
A = 2 * 3,14 / 880 = 0,007
B = 662
C = 34,1


Die nonlinea r Rückschritt Ionen konvergieren bei diesem Ansatz sehr schnell.Das Modell konvergierte in 5 Iterationen unter Verwendung dieser Werte. Siehe unten

Re: Datenanpassung an trigonometrische Funktionen FunctionSuper, genau das war die Intuition, die mir fehlte.Stattdessen habe ich die Parameterwerte in der nichtlinearen Plattform manipuliert, bis das Diagramm den Daten nahe kam, aber zu wissen, was die Parameter darstellen (und gute erste Schätzungen aus den Daten zu erhalten) ist ein viel besserer Weg.
Re: Datenanpassung an trigonometrische FunktionenHallo, vielen Dank für die Zeit, die Sie dafür aufwenden. Dies sieht nach einem großartigen Ansatz für meine Anwendung aus. aber wie hast du "automatisch" die Periodenlänge bekommen? das habe ich nicht mitbekommen..
Yoav
Re: Datenanpassung an trigonometrische FunktionenDu könntest a . verwenden Nulldurchgangsdetektor , was ein kleines bisschen JSL wäre, um durch Ihre Daten zu gehen und jedes Mal zu bemerken, wenn die Daten von negativ zu positiv wechseln. Das ist der Beginn eines Sinuswellenzyklus. Der Abstand zum nächsten Nulldurchgang ist die Periode.
Möglicherweise haben Sie Probleme mit verrauschten Daten, die mehrere Male erscheinen.
Der Nulldurchgang ist leichter zu erkennen als Spitzen.
Ich habe diese Idee in verwendet Panne .