| 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/
____________________________________________________________ ________
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
