mathematica包络检测数据平滑 [英] mathematica envelope detection data smoothing

问题描述

下面的Mathematica代码生成一个高度振荡的图.我只想绘制图的下包络线，但不知道如何绘制.任何建议将不胜感激.

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]

推荐答案

我不知道您希望它看起来如何，但是这里有一种蛮力的方法，对于我来说这是足够的起点，并且可以进一步调整:

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}]

发生的事情是，您的振荡函数具有一些非平凡的精细结构，我们需要很多要点来解决它.我们从收割-播种"图收集这些点，然后过滤掉局部最小值.由于结构良好，我们需要做两次.您实际需要的绘图存储在"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.

实际上，如果将PlotPoints的数量从50000增加到200000，然后从localMin中重复删除局部最大值，则可以获得更好的图.请注意，它将运行较慢，但是需要更多的内存.这是更改:

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}]

这是情节(做为GraphicsGrid[{{env}, {Show[{plot, env}]}}])

这篇关于mathematica包络检测数据平滑的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！