Fourier discrete transform [duplicate]

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?

marked as duplicate by BlacKow, Feyre, m_goldberg, Wjx, MarcoBOct 19 '16 at 2:54

• your data has to be 1-dimensional – BlacKow Oct 17 '16 at 22:29
• Here's the FFT plot for one of the data sets: ListLinePlot[ Abs[Fourier[Rest@Transpose[data[[1]]], FourierParameters -> {1, -1}]], PlotRange -> All] – bill s Oct 17 '16 at 22:43
• Thank you. Can you explain why did you use {1, -1} like FourierParameters? – John Oct 18 '16 at 3:29
• As you see in reference.wolfram.com/language/ref/Fourier.html a and b are parameters of the discrete Fourier transform. {-1,1} should be used in the case of data analysis, {1,-1} for signal processing. When you plot the fft you will find that the zero component corresponding to the y offset of your data is different. – mrz Oct 18 '16 at 8:26

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


• Hi Daniel, can you give some insight here? Right now I am getting the impression that I can't consider anything exact except in the trivial cases (such as an order n polynomial having n roots – I expect NSolve will always find precisely n of them). – Szabolcs Oct 26 '16 at 13:03