Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
1.0k views
in Technique[技术] by (71.8m points)

math - mathematica envelope detection data smoothing

The following Mathematica code generates a highly oscillatory plot. I want to plot only the lower envelope of the plot but do not know how. Any suggestions wouuld be appreciated.

tk0 = [Theta]'[t]*[Theta]'[t] - [Theta][t]*[Theta]''[t]
tk1 = [Theta]''[t]*[Theta]''[t] - [Theta]'[t]*[Theta]'''[t]
a = tk0/Sqrt[tk1]
f = Sqrt[tk1/tk0]
s =
 NDSolve[{[Theta]''[t] + [Theta][t] - 0.167 [Theta][t]^3 == 
    0.005 Cos[t - 0.5*0.00009*t^2], [Theta][0] == 0, [Theta]'[0] == 
    0}, [Theta], {t, 0, 1000}]

Plot[Evaluate  [f /. s], {t, 0, 1000}, 
 Frame -> {True, True, False, False}, 
 FrameLabel -> {"t", "Frequency"}, 
 FrameStyle -> Directive[FontSize -> 15], Axes -> False]

Mathematica graphics

See Question&Answers more detail:os

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Answer

0 votes
by (71.8m points)

I don't know how fancy you want it to look, but here is a brute force approach which would be good enough for me as a starting point, and can probably be tweaked further:

tk0 = [Theta]'[t]*[Theta]'[t] - [Theta][t]*[Theta]''[t];
tk1 = [Theta]''[t]*[Theta]''[t] - [Theta]'[t]*[Theta]'''[t];
a = tk0/Sqrt[tk1];
f = Sqrt[tk1/tk0];
s = NDSolve[{[Theta]''[t] + [Theta][t] - 0.167 [Theta][t]^3 == 
 0.005 Cos[t - 0.5*0.00009*t^2], [Theta][0] == 0, [Theta]'[0] ==
  0}, [Theta], {t, 0, 1000}];

plot = Plot[Evaluate[f /. s], {t, 0, 1000}, 
  Frame -> {True, True, False, False}, 
  FrameLabel -> {"t", "Frequency"}, 
  FrameStyle -> Directive[FontSize -> 15], Axes -> False];

Clear[ff];
Block[{t, x}, 
  With[{fn = f /. s}, ff[x_?NumericQ] = First[(fn /. t -> x)]]];


localMinPositionsC = 
  Compile[{{pts, _Real, 1}},
    Module[{result = Table[0, {Length[pts]}], i = 1, ctr = 0},
      For[i = 2, i < Length[pts], i++,
        If[pts[[i - 1]] > pts[[i]] && pts[[i + 1]] > pts[[i]],
          result[[++ctr]] = i]];
      Take[result, ctr]]];

(* Note: takes some time *)
points = Cases[
   Reap[Plot[(Sow[{t, #}]; #) &[ff[t]], {t, 0, 1000}, 
      Frame -> {True, True, False, False}, 
      FrameLabel -> {"t", "Frequency"}, 
      FrameStyle -> Directive[FontSize -> 15], Axes -> False, 
      PlotPoints -> 50000]][[2, 1]], {_Real, _Real}];

localMins = SortBy[Nest[#[[ localMinPositionsC[#[[All, 2]]]]] &, points, 2], First];

env = ListPlot[localMins, PlotStyle -> {Pink}, Joined -> True];

Show[{plot, env}]

What happens is that your oscillatory function has some non-trivial fine structure, and we need a lot of points to resolve it. We collect these points from Plot by Reap - Sow, and then filter out local minima. Because of the fine structure, we need to do it twice. The plot you actually want is stored in "env". As I said, it probably could be tweaked to get a better quality plot if needed.

Edit:

In fact, much better plot can be obtained, if we increase the number of PlotPoints from 50000 to 200000, and then repeatedly remove points of local maxima from localMin. Note that it will run slower and require more memory however. Here are the changes:

(*Note:takes some time*)
points = Cases[
Reap[Plot[(Sow[{t, #}]; #) &[ff[t]], {t, 0, 1000}, 
  Frame -> {True, True, False, False}, 
  FrameLabel -> {"t", "Frequency"}, 
  FrameStyle -> Directive[FontSize -> 15], Axes -> False, 
  PlotPoints -> 200000]][[2, 1]], {_Real, _Real}];

localMins =  SortBy[Nest[#[[localMinPositionsC[#[[All, 2]]]]] &, points, 2], First];

localMaxPositionsC =
  Compile[{{pts, _Real, 1}},
    Module[{result = Table[0, {Length[pts]}], i = 1, ctr = 0},
     For[i = 2, i < Length[pts], i++,
      If[pts[[i - 1]] < pts[[i]] && pts[[i + 1]] < pts[[i]], 
        result[[++ctr]] = i]];
      Take[result, ctr]]];

localMins1 = Nest[Delete[#, List /@ localMaxPositionsC[#[[All, 2]]]] &, localMins, 15];

env = ListPlot[localMins1, PlotStyle -> {Pink}, Joined -> True];

Show[{plot, env}]

Edit: here is the plot (done as GraphicsGrid[{{env}, {Show[{plot, env}]}}])

alt text


与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome to OStack Knowledge Sharing Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

...