Ajustement des données aux fonctions trigonométriquesBonjour,
J'ai du mal à ajuster les données mesurées à une fonction trigonométrique Ksin(ax+b)+c. c'est la nature de mon processus et je cherche un moyen automatique de l'adapter par code jsl.
J'ai écrit un code qui normalise les données afin que l'oscillation soit comprise entre -1 et 1.
A la fin j'aimerais avoir les constantes a et b (K et c aussi mais c'est beaucoup plus simple)
exemple ci-joint pour mon type de données.
BR,
Yoav
AnalytiqueTableau de donnéesAutreScriptRe: Ajustement des données aux fonctions trigonométriquesJe ne sais pas si c'est ce que vous voulez, mais vous pouvez résoudre ce problème en utilisant la plate-forme non linéaire.Vous trouverez ci-joint vos données, ainsi qu'une colonne de prédiction (utilisant les paramètres pour estimer) et un script qui résout ces paramètres.J'admettrai qu'il m'a fallu un certain temps pour obtenir des estimations initiales des paramètres qui ont permis à la solution de converger (probablement en raison de mon manque de familiarité avec la fonction sinus - cela fait longtemps que je n'ai pas travaillé avec de telles fonctions).Mais, conceptuellement, la plate-forme non linéaire accomplira ce que vous voulez.
Re: Ajustement des données aux fonctions trigonométriquesVoici comment obtenir une estimation de la une paramètre à l'aide de la fonction FFT. Les commentaires sur les secondes et les heures supposent que les données ont été collectées environ toutes les 4 secondes pendant environ une heure ; la FFT ne se soucie pas des unités. La FFT classe les fréquences dans les données de telle sorte que le premier bac après le bac de décalage DC représente la fréquence la plus basse qui peut faire un cycle complet dans les données d'entrée, le bac suivant est de 2 cycles complets, ... le bac du milieu est ndata/ 2 cycles. Les bacs après sont un miroir du premier.
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);
Une réponse forte dans le 4ème bac après le bac DC
Composante de fréquence la plus élevée à 0,00116448326055313 Hz ou 859 secondes par cycle
a_from_fft = 0,007316664113164 ;
La composante DC est liée à la moyenne des données, la c valeur.
Re: Ajustement des données aux fonctions trigonométriquesJe pense qu'il est utile d'utiliser ce que j'appelle la forme standard d'une fonction sinus : Y = K*sin(A*(x -
) + C.
Dans vos données :
K = 1/2 de la plage de valeurs
A = 2*pi / période
B = décalage vers la droite de 0 au point de départ de l'onde sinusoïdale
C = moyenne de vos valeurs Y
Les estimations de départ pour la régression non linéaire peuvent être trouvées à l'aide de ces définitions de la plate-forme de distribution et de Graph Builder :
K = 1/2 * (Max - Min) = 1/2* (47 - 17,6) = 14,7
A = 2 * 3,14 / 880 = 0,007
B = 662
C = 34,1


La non-linéa r régresser ion converge très rapidement en utilisant cette approche.Le modèle a convergé en 5 itérations en utilisant ces valeurs. Voir ci-dessous

Re: Ajustement des données aux fonctions trigonométriquesSuper, c'est exactement l'intuition qui me manquait.Au lieu de cela, j'ai manipulé les valeurs des paramètres dans la plate-forme non linéaire jusqu'à ce que le graphique soit proche des données, mais savoir ce que les paramètres représentent (et obtenir de bonnes estimations initiales à partir des données) est une méthode bien meilleure.
Re: Ajustement des données aux fonctions trigonométriquesBonjour, merci à tous pour le temps que vous y consacrez. Cela ressemble à une excellente approche pour mon application. mais- comment avez-vous obtenu "automatiquement" la durée de la période ? Je n'ai pas eu ça..
Yoav
Re: Ajustement des données aux fonctions trigonométriquesVous pourriez utiliser un détecteur de passage à zéro , ce qui serait un petit peu de JSL pour parcourir vos données et remarquer chaque fois que les données passent du négatif au positif. C'est le début d'un cycle d'onde sinusoïdale. La distance jusqu'au prochain passage à zéro est la période.
Vous pourriez avoir des problèmes avec des données bruitées qui semblent se croiser plusieurs fois.
Le passage à zéro est plus facile à détecter que les pics.
J'ai utilisé cette idée dans Problème .