Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

There is an elegant way to obtain the envelope of an oscillation function explained here http://mathematica.stackexchange.com/a/27854/15966

My question is how to modify the code below so I can fit data to the obtained envelope with a parameter. How do I introduce a parameter into the module so it can be used later for fitting when the result of the module is an interpolation function.

Lets' generate some data

a = RandomReal[{20, 30}];  
b = RandomReal[{5, 10}];
data = Table[{x, Cos[a*x]^2*Exp[-b*x]}, {x, 0, 1, 0.001}];
ListPlot[data, PlotRange -> All]

The code generates something like this:

data

Now the module below produces envelope for the given oscillating function

FunctionEnvelope[f_, {t_, a_, b_}, n_: 40] := 
 Module[{seeds, x, y, points, progress = 0, tempf, union}, 
  seeds = Rescale[Range[0, 1, 1/n] + 1/(2 n), {0, 1}, {a, b}];
  points = 
   Last[Last[
     Reap[Monitor[
       Do[Quiet@
         Check[progress++; {y, x} = FindMaximum[Abs[f], {t, x0}];
          x = t /. x; If[a <= x <= b, Sow[{x, y}]], Null], {x0, 
         seeds}], 
       ProgressIndicator[progress, {0, Length[seeds]}]]]]];
  union[] := 
   points = 
    Union[points, 
     SameTest -> (Abs[First[#1] - First[#2]]/
          Replace[Max[Abs[First[#1]], Abs[First[#2]]], 
           u_ /; u == 0 :> 1] < 10^-6 &)];
  union[];
  tempf = Interpolation[points];
  points = Quiet[Join[{{a, tempf[a]}}, points, {{b, tempf[b]}}]];
  union[];
  Interpolation[points]]

It can be used as follows:

q = 50;
f[t_] := E^(- 10*Abs[t]) * Cos[1000*t]/2 + 
  E^(-10*Abs[t])*Cos[(1000 + q)*t ]/2
g = FunctionEnvelope[f[x], {x, 0, 1}, 200];
Plot[{f[x], g[x]}, {x, 0, 1}, PlotRange -> All]

enter image description here

My question is how can I modify the module above to make it a function of the parameter q, so that the resulting interpolation function can be used to fit the data to get the optimal parameter q.

(The function f[t] is in principle can be moved to the definition of the module, I'm not going to change it.)

Thanks!

share|improve this question
    
Pretty sure the envelope of your function is just E^(-10 Abs[t]) Abs[Cos[q t/2]]. –  Rahul Narain Jun 18 at 9:30
    
Thank you for your reply. This is certainly true for the example. However, in general I'm going to have a function like this f[t_] := E^(- g1*Abs[t]) * Cos[1000*t]/2 + E^(-g2*Abs[t])*Cos[(1000 + q)*t ]/2 where g1 and g2 will have to be found to fit the initial data. I simplified it here for convenience, to understand how to insert a parameter into the definition of a module. –  Mikayel Jun 18 at 12:13

1 Answer 1

still unclear what you are asking - is this close?

 f[q_, t_] := E^(-10*Abs[t])*Cos[1000*t]/2 + E^(-10*Abs[t])*Cos[(1000 + q)*t]/2
 g[qq_?NumericQ] := g[qq] = Evaluate[Block[{x}, FunctionEnvelope[f[qq, x], {x, 0, 1}, 200]]];

 Plot[{f[50, x], g[50][x]}, {x, 0, 1}, PlotRange -> All]

same as yours

 Plot[Table[g[q][x], {q, {2, 10, 50}}], {x, 0, 1}, PlotRange -> All]

enter image description here

edit

the is s useful thing to do whenever having issues with fitting parameters:

 Manipulate[
     Plot[g[q][x], {x, 0, 1}, PlotRange -> {0, 1}], {{q, 50}, 48, 52}]

If you run this you'll see the function jumps around somewhat erratically for small changes in the parameter. This is why you get an error " the gradient is larger than the tolerance " when running NonlinearModelFit

share|improve this answer
    
Thanks for the reply. I mean. that the command NonlinearModelFit[data, g[q][x], q, x] does not go, also for your definition of g[q]. –  Mikayel Jun 19 at 8:15
    
fix g to take only numeric arguments (see edit) then it works with NonlinearModelFit (crummy fit though if you mean to fit to the same data ) –  george2079 Jun 19 at 16:35
    
I guess it's crummy, because it does not work -). –  Mikayel Jun 20 at 8:13
    
I guess it's crummy, because it does not work -). Let's try data = Table[{x, Abs[Cos[44*x/2]]*Exp[-10*x]}, {x, 0, 1, 0.001}]; This generates data, that fits the envelope with q=44 very well (check Show[Plot[{g[44][x]}, {x, 0, 1}, PlotRange -> All, PlotStyle -> {Thick, Red}], ListPlot[data, PlotRange -> All]] ). If we now let NonlinearModelFit try to find it near q=40: nlm = NonlinearModelFit[data, g[q][x], {{q, 40}}, x] nlm["BestFitParameters"] with your new definition of g[q], it does not work: NonlinearModelFit::sszero: The step size in the search has become.. {q -> 41.} –  Mikayel Jun 20 at 8:21
    
it works in the sense that it is in a form suitable for input to NonlinearSolve. The fact you don't get a good result has to do with the behavior of your function, which is an entirely different issue. There are a bunch of questions on this site dealing with fitting difficulties. Typically you need to get very close with an initial guess at the solution. –  george2079 Jun 20 at 12:02

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.