Adattamento dei dati alle funzioni trigonometricheCiao,
Sto lottando con l'adattamento dei dati misurati a una funzione trigonometrica Ksin (ax + b) + c. questa è la natura del mio processo e sto cercando un modo automatico per adattarlo al codice jsl.
Ho scritto un codice che normalizza i dati in modo che l'oscillazione sia tra -1 e 1.
Alla fine vorrei ottenere le costanti a e b (anche K e c ma è molto più semplice)
esempio allegato per il mio tipo di dati.
BR,
Yoav
AnalisiTabella datiAltroScriptingRe: Adattamento dei dati alle funzioni trigonometricheNon so se questo è quello che vuoi, ma puoi risolverlo usando la piattaforma non lineare.In allegato ci sono i tuoi dati, insieme a una colonna predittore (usando i parametri per stimare) e uno script che risolve per quei parametri.Devo ammettere che mi ci è voluto un po' di tempo per ottenere ipotesi iniziali per i parametri che hanno permesso alla soluzione di convergere (probabilmente a causa della mia mancanza di familiarità con la funzione seno - è passato molto tempo da quando ho lavorato con tali funzioni).Ma, concettualmente, la piattaforma non lineare realizzerà ciò che desideri.
Re: Adattamento dei dati alle funzioni trigonometricheEcco come ottenere una stima del un parametro utilizzando la funzione FFT. I commenti su secondi e ore presumono che i dati siano stati raccolti circa ogni 4 secondi per circa un'ora; alla FFT non interessano le unità. La FFT raggruppa le frequenze nei dati in modo tale che il primo bin dopo il bin di offset DC rappresenti la frequenza più bassa che può fare un ciclo completo nei dati di input, il bin successivo è 2 cicli completi, ... il bin centrale è ndata/ 2 cicli. I cassonetti dopo sono uno specchio del primo.
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);
Una risposta forte nel quarto bin dopo il bin DC
Componente di frequenza più grande a 0,00116448326055313 Hz o 859 secondi per ciclo
a_from_fft = 0.007316664113164;
La componente DC è correlata alla media dei dati, il C valore.
Re: Adattamento dei dati alle funzioni trigonometrichePenso che sia utile usare quella che chiamo la forma standard di una funzione seno: Y = K*sin(A*(x -
) + C.
Nei tuoi dati:
K = 1/2 dell'intervallo di valori
A = 2*pi greco / punto
B = spostamento a destra da 0 al punto iniziale dell'onda sinusoidale
C = media dei tuoi valori Y
Le stime di partenza per la regressione non lineare possono essere trovate utilizzando queste definizioni dalla piattaforma di distribuzione e 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 nonlinea regredire ione converge molto rapidamente utilizzando questo approccio.Il modello converge in 5 iterazioni utilizzando questi valori. Vedi sotto

Re: Adattamento dei dati alle funzioni trigonometricheOttimo, è proprio l'intuizione che mi mancava.Invece, ho manipolato i valori dei parametri nella piattaforma non lineare fino a quando il grafico non sembrava vicino ai dati, ma sapere cosa rappresentano i parametri (e ottenere buone ipotesi iniziali dai dati) è un modo molto migliore per andare.
Re: Adattamento dei dati alle funzioni trigonometricheCiao, grazie a tutti per il tempo che ci dedichi. Questo sembra un ottimo approccio per la mia applicazione. ma- come hai ottenuto "automaticamente" la durata del periodo? non l'ho capita..
Yoav
Re: Adattamento dei dati alle funzioni trigonometrichePotresti usare un rilevatore di passaggio per lo zero , che sarebbe un po' di JSL per esaminare i tuoi dati e notare ogni volta che i dati passano da negativi a positivi. Questo è l'inizio di un ciclo di onde sinusoidali. La distanza dal prossimo passaggio per lo zero è il periodo.
Potresti avere problemi con dati rumorosi che sembrano incrociarsi più volte.
Il passaggio per lo zero è più facile da rilevare rispetto ai picchi.
Ho usato quell'idea in Glitch .