сглаживание данных обнаружения конвертов Mathematica

Question

Следующий код Mathematica генерирует сильно колебательный график. Я хочу построить только нижнюю оболочку сюжета, но не знаю как. Будем благодарны за любые предложения.

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]

Answer 1

Я не знаю, насколько вы хотите, чтобы это выглядело, но вот подход грубой силы, который был бы достаточно хорош для меня в качестве отправной точки, и, вероятно, может быть изменен дальше:

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

Что происходит, так это то, что ваша колебательная функция имеет некоторую нетривиальную тонкую структуру, и нам нужно много точек для ее решения. Мы собираем эти точки из Plot by Reap - Sow, а затем отфильтровываем локальные минимумы. Из-за тонкой структуры нам нужно сделать это дважды. Сюжет, который вы на самом деле хотите, хранится в «env». Как я уже сказал, возможно, его можно настроить, чтобы получить более качественный сюжет, если это необходимо.

Edit:

Фактически, намного лучший график может быть получен, если мы увеличим количество PlotPoints с 50000 до 200000, а затем многократно удаляем точки локальных максимумов из localMin. Обратите внимание, что он будет работать медленнее и, тем не менее, потребовать больше памяти. Вот изменения:

(*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}]}}])

Answer 2

Решение на основе изображений

Я не претендую на это, ни крепкий, ни общий. Но это быстро и весело. Он использует преобразования изображения, чтобы найти края (возможно, из-за сильного колебательного характера вашей функции):

Функция:

envelope[plot_] := Module[{boundary, Pr, rescaled},

  (* "rasterize" the plot, identify the lower edge and isolate pixels*)

  boundary = Transpose@ImageData@Binarize@plot /. {x___, 0, 1, y___} :>
     Join[Array[1 &, Length[{x}]], {0}, Array[1 &, Length[{y}] + 1]];

  (* and now rescale *)

  Pr = PlotRange /. Options[plot, PlotRange];
  rescaled = Position[boundary, 0] /.
    {x_, y_} :> {
      Rescale[x, {1, Dimensions[boundary][[1]]}, Pr[[1]]],
      Rescale[y, {1, Dimensions[boundary][[2]]}, Reverse[Pr[[2]]]]
      };

  (* Finally, return a rescaled and slightly smoothed plot *)

  Return[ListLinePlot@
    Transpose@{( Transpose[rescaled][[1]])[[1 ;; -2]], 
      MovingAverage[Transpose[rescaled][[2]], 2]}]
   ]  

Код тестирования:

tk0 = phi'[t] phi'[t] - phi[t] phi''[t];
tk1 = phi''[t] phi''[t] - phi'[t] phi'''[t];
a = tk0/Sqrt[tk1];
f = Sqrt[tk1/tk0];
s = NDSolve[{
    phi''[t] + phi[t] - 0.167 phi[t]^3 == 
     0.005 Cos[t - 0.5*0.00009*t^2],
    phi[0] == 0,
    phi'[0] == 0},
   phi, {t, 0, 1000}];
plot = Plot[Evaluate[f /. s], {t, 0, 1000}, Axes -> False];
Show[envelope[plot]]

Редактировать

Исправляя ошибку в коде выше, результаты более точны:

envelope[plot_] := Module[{boundary, Pr, rescaled},

  (*"rasterize" the plot,
  identify the lower edge and isolate pixels*)

  boundary = 
   Transpose@ImageData@Binarize@plot /. {x___, 0, 1, y___} :> 
     Join[Array[1 &, Length[{x}]], {0}, Array[1 &, Length[{y}] + 1]];

  (*and now rescale*)

  Pr = PlotRange /. Options[plot, PlotRange];

  rescaled = Position[boundary, 0] /. {x_, y_} :>
     {Rescale[
       x, {(Min /@ Transpose@Position[boundary, 0])[[1]], (Max /@ 
           Transpose@Position[boundary, 0])[[1]]}, Pr[[1]]], 
      Rescale[y, {(Min /@ 
           Transpose@Position[boundary, 0])[[2]], (Max /@ 
           Transpose@Position[boundary, 0])[[2]]}, Reverse[Pr[[2]]]]};

  (*Finally,return a rescaled and slightly smoothed plot*)
  Return[ListLinePlot[
    Transpose@{(Transpose[rescaled][[1]])[[1 ;; -2]], 
      MovingAverage[Transpose[rescaled][[2]], 2]}, 
    PlotStyle -> {Thickness[0.01]}]]]

, .