| mpfit of parametric data? [message #40290] |
Wed, 28 July 2004 08:45  |
jamiesmyth_uni@yahoo.
Messages: 6
Registered: July 2004
|
Junior Member |
|
|
I suppose this is really a question for Craig but I figure here is as
good a place to ask as any... Does anyone know how I can go about
fitting parametric data using MPFIT? I have done a fair bit of 1d
fitting with mpfit (MPFITEXPR) but I'm really stumped on this one. I
want to fit to the following parametric parameterisation:
x = (a0+a1*t) + sin(a2*t+phase)
y = (b0+b1*t) + cos(b2*t+phase)
where, a0, a1, a3, b0, b1, b2 and phase are all fit parameters.
The intention is to try and fit the motion of a spinning top that
precesses. I have very long running (but noisy) time series data for
the x and y values. Alternatively, you can think of me having x(t) and
y(t) sampled at identical times. I am mainly interested in the phase
parameter.
This is proving considerably more difficult than I expected it to be!
Thanks.
Jamie
|
|
|
|
| Re: mpfit of parametric data? [message #40508 is a reply to message #40290] |
Mon, 16 August 2004 11:14  |
Craig Markwardt
Messages: 1869
Registered: November 1996
|
Senior Member |
|
|
"jamiesmyth_uni@yahoo.ca" <jamiesmyth_uni@yahoo.ca> writes:
> After two weeks of vacation I'm back at this. I've de-trended the
> original time series by fitting to a quadratic, and estimated the
> frequencies of the components by looking at the power spectrum.
> Unfortunately, I still cannot fit the amplitude and phase of a trivial
> sinusoid such as 'A*sin(2*!dpi*w*t+phi)'. How do I go about estimating
> the phase of the following trivial example?
>
> IDL> n = 2048
> IDL> t = dindgen(n) * 0.125 ; time
> IDL> freq = dindgen(n)/(n*0.128)
> IDL> p0 = [0.03, 0.06923, 2.3]
> IDL> data = p0(0)*sin( 2*!dpi*p0(1)*t + p0(2) )
> IDL> plot, t, data
> IDL> Ft_data = fft(data)
> IDL> plot, freq, abs(Ft_data)^2, xrange=[0,0.5] ; frequency estimate
> IDL> plot, freq, atan(double(Ft_data),imaginary(Ft_data)),
> xrange=[0,0.5] ; ?phase estimate?
>
> I think I understand what Craig said about a local-minimum but I'm a
> little surprised that such a simple problem (i.e. estimating the
> amplitude, frequency and phase of a sine wave) would be so difficult?
> How is it that my Lecroy/Tektronix scope can solve this in real time
> but I cannot do it with IDL and a dual xeon? Surely I must be missing
> something...
Jamie, since you don't have any call to MPFIT here, I am a little
confused about how your question relates to MPFIT.
Also, you have a mismatch in your FREQ and T variables, right?
The 2-argument version of ATAN accepts variables as ATAN(Y,X), not
ATAN(X,Y), so shouldn't you swap the order of the imaginary and real
components?
Finally, you are plotting the spectrum of phases. Don't you want the
phase at the maximum Fourier power? And, how do you know the phase is
wrong?
Craig
--
------------------------------------------------------------ --------------
Craig B. Markwardt, Ph.D. EMAIL: craigmnet@REMOVEcow.physics.wisc.edu
Astrophysics, IDL, Finance, Derivatives | Remove "net" for better response
------------------------------------------------------------ --------------
|
|
|
|
| Re: mpfit of parametric data? [message #40512 is a reply to message #40290] |
Mon, 16 August 2004 08:36 |
Paolo Grigis
Messages: 171
Registered: December 2003
|
Senior Member |
|
|
When you get your starting guess for the frequency (let's call it
guessfreq) from the FFT you should be able to use mpfitfun with
starting parameters, say, [1.,guessfreq,0.] and obtain a good fit...
IDL> dummy=max(abs(Ft_data),count)
IDL> guessfreq=freq[count]
IDL> print,guessfreq
0.068664548
IDL> par=mpfitfun('yoursinusoid',t,data,errors,[1.,guessfreq,0.])
Iter 1 CHI-SQUARE = 1.0754572E+09
P(0) = 1.00000
P(1) = 0.0686645
P(2) = 0.00000
[...]
Iter 7 CHI-SQUARE = 1.1268491E-23
P(0) = -0.0300000
P(1) = 0.0692300
P(2) = -0.841593
You may obtain a negative amplitude and/or phase sometimes, but if
this bothers you, you can use the PARINFO structure in mpfitfun for
constraining them to positive values.
In this simple case it looks like you can fit your data
without needing good guesses on the amplitude or the phase.
But possibly more work is needed for the real data...
Ciao,
Paolo
jamiesmyth_uni@yahoo.ca wrote:
> After two weeks of vacation I'm back at this. I've de-trended the
> original time series by fitting to a quadratic, and estimated the
> frequencies of the components by looking at the power spectrum.
> Unfortunately, I still cannot fit the amplitude and phase of a trivial
> sinusoid such as 'A*sin(2*!dpi*w*t+phi)'. How do I go about estimating
> the phase of the following trivial example?
>
> IDL> n = 2048
> IDL> t = dindgen(n) * 0.125 ; time
> IDL> freq = dindgen(n)/(n*0.128)
> IDL> p0 = [0.03, 0.06923, 2.3]
> IDL> data = p0(0)*sin( 2*!dpi*p0(1)*t + p0(2) )
> IDL> plot, t, data
> IDL> Ft_data = fft(data)
> IDL> plot, freq, abs(Ft_data)^2, xrange=[0,0.5] ; frequency estimate
> IDL> plot, freq, atan(double(Ft_data),imaginary(Ft_data)),
> xrange=[0,0.5] ; ?phase estimate?
>
> I think I understand what Craig said about a local-minimum but I'm a
> little surprised that such a simple problem (i.e. estimating the
> amplitude, frequency and phase of a sine wave) would be so difficult?
> How is it that my Lecroy/Tektronix scope can solve this in real time
> but I cannot do it with IDL and a dual xeon? Surely I must be missing
> something...
>
> Any help is greatly appreciated!
> Jamie
>
--
____________________________________________________________ ________
Paolo Grigis
ETHZ - Institute of Astronomy email: pgrigis@astro.phys.ethz.ch
Scheuchzerstrasse 7
ETH Zentrum phone: ++41 1 632 42 20
8092 Zurich fax : ++41 1 632 12 05
Switzerland http://www.astro.phys.ethz.ch/
____________________________________________________________ ________
|
|
|
|
| Re: mpfit of parametric data? [message #40533 is a reply to message #40290] |
Thu, 12 August 2004 15:00 |
jamiesmyth_uni@yahoo.
Messages: 6
Registered: July 2004
|
Junior Member |
|
|
After two weeks of vacation I'm back at this. I've de-trended the
original time series by fitting to a quadratic, and estimated the
frequencies of the components by looking at the power spectrum.
Unfortunately, I still cannot fit the amplitude and phase of a trivial
sinusoid such as 'A*sin(2*!dpi*w*t+phi)'. How do I go about estimating
the phase of the following trivial example?
IDL> n = 2048
IDL> t = dindgen(n) * 0.125 ; time
IDL> freq = dindgen(n)/(n*0.128)
IDL> p0 = [0.03, 0.06923, 2.3]
IDL> data = p0(0)*sin( 2*!dpi*p0(1)*t + p0(2) )
IDL> plot, t, data
IDL> Ft_data = fft(data)
IDL> plot, freq, abs(Ft_data)^2, xrange=[0,0.5] ; frequency estimate
IDL> plot, freq, atan(double(Ft_data),imaginary(Ft_data)),
xrange=[0,0.5] ; ?phase estimate?
I think I understand what Craig said about a local-minimum but I'm a
little surprised that such a simple problem (i.e. estimating the
amplitude, frequency and phase of a sine wave) would be so difficult?
How is it that my Lecroy/Tektronix scope can solve this in real time
but I cannot do it with IDL and a dual xeon? Surely I must be missing
something...
Any help is greatly appreciated!
Jamie
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
