Elegant way of obtaining the envelope of oscillating function [duplicate]

Question

I am solving a differential equation numerically and the output is an oscillating function with the amplitude of the oscillation decaying in time. I would like to extract the power law governing this amplitude.

Let's make this concrete. Let's say I am solving the system

x''[t] + 2t^(-1) x'[t] + t^(-2) x[t] == 0

which has solutions

x= C[1] t^(-1/2) Cos[Sqrt[3]/2 Log[t]]

and the corresponding $\sin$. I would like to extract the power law $t^{-1/2}$.

Of course, my system is much more complicated and I have to solve it numerically. I get as my solution an interpolating function. I was thinking of generating a table of function values, sorting according to the maxima and then fitting a power law, but maybe there is a better way?

Thanks!

Answer 1

Here is a general envelope method developed for this SystemModeler industry example:

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's used like this:

f[x_] := x^2 Sin[x] Sin[x^2]

g = FunctionEnvelope[f[x], {x, 0, 15}, 100];

Now g[x] is the automatically determined envelope function:

Plot[{f[x], g[x], -g[x]}, {x, 0, 15}]

This works by smoothly interpolating between the local maxima of the function.

Answer 2

The idea of an envelope is often quite clear, but it is a slippery thing to try and define it mathematically. In audio, the envelope is often used to help locate start and ending points for events, and is often associated with two parameters: a rise time (how fast the envelope is allowed to grow) and a decay time (how fast the envelope is allowed to die away. Concretely, one can calculate the envelope as a combination of two linear systems with two different time constants.

decay = 0.006; rise = 0.2;
filt[z_, u_] := Max[decay z + (1 - decay) u, rise z + (1 - rise) u];
env = FoldList[filt, 0, Abs[data]];

The filt function implements the two linear systems and combines them with Max. The data is a list (in audio it might be a single channel extracted from a wav file), and the two parameters (both must be between zero and one) specify the rise and decay times.

Here is this enveloping operation in action. The light red is the data and the dark blue is the calculated envelope.