Fourier discrete transform

Question

Code is given below

m[t_] := {mx[t], my[t], mz[t]}

γ = 28;
h = 6.62*10^-34;
e = 1.6*10^-19;

Subscript[μ, 0] = 1.25*10^-6;
Subscript[μM, 0] = 800*10^-3;
Subscript[M, 0] = 0.64*10^6;
Subscript[r, 0] = 100*10^-9;
Subscript[l, 0] = 3*10^-9;
Subscript[I, dc] = 1*10^-3;
Subscript[B, dc] = 200*10^-3;
Subscript[α, G] = 0.01;

p = {0, 0, 1};
σ =(γ*h/2*e)*1/(Subscript[M, 0]*Pi*(Subscript[r, 0])^2)*Subscript[l, 0];
Subscript[B, eff] = {Subscript[B, dc], 0, 0}-Subscript[μM, 0]*(m[t]*p);

system1 ={D[m[t], t] ==γ*(Cross[Subscript[B, eff], m[t]]) + Subscript[α, 
G]*(Cross[m[t], D[m[t], t]]) +σ*Subscript[I, dc]*(Cross[m[t], Cross[m[t], 
p]]),(m[t] /. t -> 0) == {0, 1, 0}};

s1 = NDSolve[system1, m[t], {t, 0, 50}]

Plot[Evaluate[{mx[t], my[t], mz[t]} /. s1], {t, 0, 50},AxesLabel -> {t, m}]
Plot[Evaluate[mx[t] /. s1], {t, 0, 5}, AxesLabel -> {t, mx}]
z = Reap[s1 = NDSolve[{system1, WhenEvent[mx'[t] == 0, Sow[{t, mx[t]}]]}, 
m[t], {t, 0, 5}]][[2, 1]]
data = Table[z, {t, 0, 5}];
ListLinePlot[data, Joined -> True]
ListLinePlot[Abs[Fourier[data]], Joined -> True, PlotRange -> Full]

I need to find the frequency using the Fourier discrete transform. So first I need to do this transform, but I got a strange result. What am I doing wrong?

Answer 1

With everything set up as in the question, one can estimate a frequency just by counting peaks in some range and dividing. The code below just slightly changes code from the post.

z = Reap[s1 = 
     NDSolve[{system1, WhenEvent[mx'[t] == 0, Sow[{t, mx[t]}]]}, 
      m[t], {t, 0, 3}]][[2, 1]];
Pi*(Length[z] - 1)/(z[[-1, 1]] - z[[1, 1]])

(* Out[373]= 21.8156118886 *)

If we go further out, say to 4 or 5, the estimate goes up. Presumably the peaks start to get closer. We can actually see this as below. I adapt some prior MSE code to show what Prony's method can do to estimate the initial part of this function. I keep singular values that are at least 1/10 of the largest. More than that and the computation seems to become numerically unstable (though maybe there are better ways to do it).

delta = .004;
xmax = 4;
vals = Table[{t, mx[t] /. s1[[1]]}, {t, 0, xmax, delta}];
nn = Length[vals];
mat = Most[Partition[vals[[All, 2]], Floor[nn/2], 1]];
keep = Length[SingularValueList[mat, Tolerance -> 10^(-1)]]
mat2 = Most[Partition[vals[[All, 2]], keep, 1]];
rhs = Drop[vals[[All, 2]], keep];
soln = PseudoInverse[mat2].rhs;
roots = xx /. NSolve[xx^keep - soln.xx^Range[0, keep - 1] == 0, xx]
freqs = Log[roots]

(* Out[420]= 3

Out[424]= {0.995719653834 - 0.0875079462907 I, 
 0.995719653834 + 0.0875079462907 I, 1.00072750797}

Out[425]= {-0.000442560954337 - 
  0.0876589023872 I, -0.000442560954337 + 
  0.0876589023872 I, 0.000727243468077} *)

newmat = Map[roots^# &, Range[0, nn - 1]];
coeffs = PseudoInverse[newmat].vals[[All, 2]];
newf = Chop[TrigExpand[ExpToTrig[coeffs.Exp[freqs*t/sr]]]]

 (* Out[428]= -0.203137465804 Cos[21.9147255968 t] Cosh[
   0.110640238584 t] + 0.404894655937 Cosh[0.181810867019 t] - 
 0.200611197705 Cosh[0.110640238584 t] Sin[21.9147255968 t] + 
 0.203137465804 Cos[21.9147255968 t] Sinh[0.110640238584 t] + 
 0.200611197705 Sin[21.9147255968 t] Sinh[0.110640238584 t] + 
 0.404894655937 Sinh[0.181810867019 t] *)

Notice that the sinusoidal parts are coming in at close to 22. A plot will show also that eventually the actual signal starts to oscillate a bit faster than this approximate reconstruction.

p1 = Plot[newf, {t, 0, xmax}, PlotStyle -> {Red}];
p2 = ListPlot[vals];
Show[{p2, p1}]