4
$\begingroup$

I'm solving a vibration problem. After solving two systems of ODEs (one valid in 0-0.1 and the other in 0.1-1) I get the plot of 3 functions, one of them is this, the plot of the variable $y_w$ (Image 1).

enter image description here

I want to see the natural frequencies of this vibration using Fourier transform and I used this code

{

setyw = {};
step = 0.0001;

While[ti < tf, ti = ti + step; 
  If[ti >= 0.1, 
   AppendTo[setyw, Evaluate[(yw[ti] /. sfreeNL)[[1, 1]]]], 
   AppendTo[setyw, Evaluate[(yw[ti] /. sNL)[[$1$]]]]]];

ListLinePlot[Abs[Fourier[setyw]], PlotRange -> Full]}

What I get from 0 to 1 with 0.0001 step is image 2. Which is totally not what I would expect, especially because if I just take a shorter interval, like from 0 to 0.5 I get image 3, which is exactly the same as before, one peak at the beginning and one peak at the end of the interval.

If I add some noise with this code

{

coeff = 1;

mag = 2;

setyw = {};

step = 0.01;

While[ti < tf, ti = ti + step; 
  If[ti >= 0.1, 
   AppendTo[setyw, 
    Evaluate[(yw[ti] /. sfreeNL)[[1, 1]]] + 
     coeff 10^-mag (RandomReal[] - 1/2)], 
   AppendTo[setyw, 
    Evaluate[(yw[ti] /. sNL)[[$1$]]] + 
     coeff 10^-mag (RandomReal[] - 1/2)]]];
ListLinePlot[Abs[Fourier[setyw]], PlotRange -> Full]}

I get image 4 and 5, which as you can see are symmetric in respect to the half of the interval.

How would you explain this? Is it a problem of code or maths?


Complete code:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
Clear["Global`*"];

mw = 40 ;

mv = 20; mp = 12; kw11 = 300000; kw12 = 300000; kw31 = 2000;

cw11 = 300; cw12 = 300; cw31 = 200;

lp = 1.3; lv = 1.4; L = 0.5; A = 0.3; H = 0.45; rp = 0.15; w = 0.01;

rv = 0.16; xG = (L^2 + A^2 + L A)/(3 (L + A)); 

yG = (H (2 L + A))/(3 (L + A));

lw2 = xG; lw1 = L + A - xG; hw = yG; hv = H - yG; P = 4000;

g = 9.81; M = rp P;

Iv = mv/12 (3 (rv^2 + (rv - w)^2) + lv^2);

Ipx = (mp rp^2)/2;
Ipz = mp/12 (3 rp^2 + lp^2);
Iw = mw (hw + hv)/(12) ((hw + hv)^2 (4 L + A) + 4 L^3 + 6 L^2 A + 4 L A^2 +A^3);

static = Solve[{ywst (kw11 lw1 - 
     kw12 lw2) + ϑwst (kw11 lw1^2 + kw12 lw2^2 - 
     kw31) - kw31 ϑvst - (mv + 
     mp) g hv ϑwst == 0,
ywst (kw11 + kw12) + ϑwst (kw11 lw1 - 
     kw12 lw2) + (mw + mv + mp) g == 0,
kw31 ϑwst + kw31 ϑvst == 
 0}, {ywst, ϑwst, ϑvst}];

sNL = NDSolve[{(mw + mv + mp) yw''[
   t] == (mv + mp) hv ϑw''[t] Sin[ϑw[t]] -
   mp ξp[t] ϑw''[
    t] Cos[ϑv[t]] + (mv + mp) hv (ϑw'[
     t])^2 Cos[ϑw[t]] + 
  mp ξp[t] (ϑv'[t])^2 Sin[ϑv[t]] - 
  mp ξp'[t] ϑv'[
    t] Cos[ϑv[
     t]] - (yw[t] + (ywst /. static)[[1]]) (kw11 + 
     kw12) - (ϑw[t] + (ϑwst /. static)[[
      1]]) (lw1 kw11 - lw2 kw12) + 
  yw'[t] (cw11 + cw12) + ϑw'[
    t] (cw11 lw1 - cw12 lw2) - (mw + mv + mp) g,
(Iw + Ipz + Iv + (mp + mv) hv^2) ϑw''[
   t] == -(Ipz + 
      mp hv ξp[
        t] (Sin[ϑv[t]] Cos[ϑw[t]] - 
         Sin[ϑw[t]] Cos[ϑv[
            t]])) ϑv''[t] + (mv + mp) hv yw''[
    t] Sin[ϑw[t]] + 
  mp hv ξp''[
    t] (Cos[ϑv[t]] Cos[ϑw[t]] + 
     Sin[ϑw[t]] Sin[ϑv[t]]) - 
  2 mp hv ξp'[
    t] (ϑw'[
       t] Sin[ϑw[t]] Cos[ϑv[
        t]] - ϑv'[
       t] Sin[ϑw[t]] Cos[ϑv[t]]) - 
  mp hv ξp[
    t] (ϑv'[
     t])^2 (Cos[ϑv[t]] Cos[ϑw[t]] - 
     Sin[ϑw[t]] Sin[ϑv[t]]) - (yw[
      t] + (ywst /. static)[[1]]) (kw11 lw1 - 
     kw12 lw2) - (ϑw[
      t] + (ϑwst /. static)[[1]]) (lw1^2 kw11 + 
     lw2^2 kw12 + kw31) + 
  kw31 (ϑv[t] + (ϑvst /. static)[[1]]) - 
  yw'[t] (-cw11 lw1 + cw12 lw2) + ϑw'[
    t] (cw11 lw1^2 + cw12 lw2^2 + cw31) - 
  cw31 ϑv'[
    t] + (mv + 
     mp) g hv Sin[ϑw[
      t] + (ϑwst /. static)[[1]]],
(Iv + Ipz + mp (ξp[t])^2) ϑv''[
   t] == -mp ξp[t] yw''[
    t] Cos[ϑv[t]] - (Iv + Ipz + 
     mp hv ξp[
       t] (Cos[ϑw[t]] Sin[ϑv[t]] - 

        Sin[ϑw[t]] Cos[ϑv[
           t]])) ϑw''[t] - 
  2 mp ξp[t] ξp'[t] ϑv'[t] - 
  mp hv ξp[
    t] (ϑw'[
     t])^2 (Cos[ϑv[t]] Cos[ϑw[t]] + 
     Sin[ϑw[t]] Sin[ϑv[t]]) - 
  kw31 (ϑv[t] + (ϑvst /. static)[[
      1]] - (ϑwst /. static)[[1]] - ϑw[
      t]) + cw31 (ϑv'[t] - ϑw'[t]) - 
  mp g ξp[
    t] Cos[ϑv[t] + (ϑvst /. static)[[
      1]]],
mp ξp''[t] == 
 mp hv ϑv''[
    t] (Cos[ϑv[t]] Cos[ϑw[t]] + 
     Sin[ϑw[t]] Sin[ϑv[t]]) - 
  mp yw''[t] Sin[ϑw[t]] - 
  mp hv (ϑw'[
     t])^2 (Sin[ϑw[t]] Cos[ϑv[t]] - 
     Sin[ϑv[t]] Cos[ϑw[t]]) - 
  mp ξp[t] (ϑv'[t])^2 - 
  mp g Sin[ϑv[t] + (ϑvst /. static)[[
      1]]] + P, yw[0] == (ywst /. static)[[1]], 
yw'[0] == 
 0, ϑw[0] == (ϑwst /. static)[[
  1]], ϑw'[0] == 0, ϑv[0] == 
 0.1, ϑv'[0] == 0, ξp[0] == 0, ξp'[0] == 
 0}, {yw, ϑw, ϑv, ξp}, {t, 0, 0.1}, 
Method -> {"EquationSimplification" -> "Residual"}, MaxSteps -> Infinity];

yw1 = Evaluate[yw[0.1] /. sNL];
[CurlyTheta]w1 = Evaluate[ϑw[0.1] /. sNL];
[CurlyTheta]v1 = Evaluate[ϑv[0.1] /. sNL];
dyw1 = Evaluate[yw'[0.1] /. sNL];
dϑw1 = Evaluate[ϑw'[0.1] /. sNL];
dϑv1 = Evaluate[ϑv'[0.1] /. sNL];

sfreeNL = 
NDSolve[{(mw + mv) yw''[t] + (mw + mv) (lw1 - lw2) ϑw''[
    t] - (mw + mv) g + (kw11 + kw12) yw[
    t] + (lw1 kw11 - lw2 kw12) ϑw[
    t] + (cw11 + cw12) yw'[
    t] + (lw1 cw11 - lw2 cw12) ϑw'[t] == 0,
(mw + mv) (lw1 - lw2) yw''[
    t] + ((mw + mv) (lw1 - lw2)^2 + Iw + Iv + 
     mv hv^2) ϑw''[
    t] - (mw + mv) g (lw1 - lw2) + (kw11 lw1 - kw12 lw2) yw[
    t] + (kw11 lw1^2 + kw12 lw2^2) ϑw[t] - 
  kw31 (ϑv[t] - ϑw[t]) + (cw11 lw1 - 
     cw12 lw2) yw'[
    t] + (cw11 lw1^2 + cw12 lw2^2 + cw31) ϑw'[t] - 
  cw31 ϑv'[t] == 0,
Iv ϑv''[t] + 
  kw31 (ϑv[t] - ϑw[t]) + 
  cw31 (ϑv'[t] - ϑw'[t]) == 0,
yw[0.1] == yw1, 
yw'[0.1] == 
 dyw1, ϑw[0.1] == ϑw1, ϑw'[
  0.1] == dϑw1, ϑv[
  0.1] == ϑv1, ϑv'[0.1] == 
 dϑv1}, {yw, ϑw, ϑv}, {t, 0.1,
 1}, Method -> {"ExplicitRungeKutta"}, MaxSteps -> Infinity];

p1 = Plot[Evaluate[yw[x] /. sNL], {x, 0, 0.1}, PlotRange -> All, 
AxesLabel -> {t, Subscript[y, w]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

p2 = Plot[Evaluate[ϑw[x] /. sNL], {x, 0, 0.1}, 
PlotRange -> All, AxesLabel -> {t, Subscript[ϑ, w]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

p3 = Plot[Evaluate[ϑv[x] /. sNL], {x, 0, 0.1}, 
PlotRange -> All, AxesLabel -> {t, Subscript[ϑ, v]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

p4 = Plot[Evaluate[ξp[x] /. sNL], {x, 0, 0.1}, PlotRange -> All, 
AxesLabel -> {t, Subscript[ξ, p]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

p12 = Plot[Evaluate[yw[x] /. sfreeNL], {x, 0.1, 1}, PlotRange -> All, 
AxesLabel -> {t, Subscript[y, w]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

p22 = Plot[Evaluate[ϑw[x] /. sfreeNL], {x, 0.1, 1}, 
PlotRange -> All, AxesLabel -> {t, Subscript[ϑ, w]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

p32 = Plot[Evaluate[ϑv[x] /. sfreeNL], {x, 0.1, 1}, 
PlotRange -> All, AxesLabel -> {t, Subscript[ϑ, v]}, 
LabelStyle -> Directive[Red, Medium], PlotStyle -> Thick];

Show[p1, p12]
Show[p2, p22]
Show[p3, p32]

ti = Input["starting time? (minimo 0)"];

tf = Input["final time? (max 1)"];

setyw = {};

step = Input["which step?"];

While[ti < tf, ti = ti + step; 
If[ti >= 0.1, 
AppendTo[setyw, Evaluate[(yw[ti] /. sfreeNL)[[1, 1]]]], 
AppendTo[setyw, Evaluate[(yw[ti] /. sNL)[[1]]]]]];

ListLinePlot[Abs[Fourier[setyw]], PlotRange -> Full]

Here is a link to the file:

https://drive.google.com/file/d/0B-w5qQPBHkwLZEFyYTU4blZTN3M/view?usp=sharing

$\endgroup$

marked as duplicate by Jens, C. E., m_goldberg, Yves Klett, Young Aug 15 '16 at 0:34

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • 1
    $\begingroup$ You can select code and press the {} button to format it. I think posting the code here without input dialogues is handier. $\endgroup$ – Feyre Aug 14 '16 at 12:15
  • $\begingroup$ @Feyre I think it's ok now isn't it? $\endgroup$ – mtlvc0 Aug 14 '16 at 12:45
  • $\begingroup$ I think you have accidentally deleted some relevant code, it doesn't run. $\endgroup$ – Feyre Aug 14 '16 at 12:47
  • $\begingroup$ @Feyre I think now is ok. I Uploaded code and a link to the file. $\endgroup$ – mtlvc0 Aug 14 '16 at 14:59
  • 2
    $\begingroup$ This is the expected output from Fourier. The negative frequencies come after the positive. This is fine because the numerical Fourier transform is periodic. I put some notes together on this here $\endgroup$ – Hugh Aug 14 '16 at 22:49
6
$\begingroup$

First, having a dataset of length $N$ sampled with a time step $\Delta t$, the Fourier (FFT) can resolve frequencies $\omega_k$ from the interval $(-\omega_{Nyq},\omega_{Nyq})$, where $\omega_{Nyq}=\frac{1}{2\Delta t}$ is the Nyquist frequency, and the frequencies are $\omega_k=\frac{k}{T}$, where $T$ is the length of the time series and $k$ is an integer from the range $(-N/2,N/2)$.

EDIT: When you apply FFT to finite data, you'll always achieve an uncomparable power of the frequency related to the length of the time series - you should simply ignore it. Judging from your image 1, there's quite a lot of periodicity, and the length of the time series is very large. In such a way you are checking also a huge number of frequencies that are not physically present in the data - their powers will be nearly infinitesimal due to finiteness of the data and discreteness of the sampling. You are really interested in a very small portion of the spectrum starting from a bit more than zero (i.e., 1 or 2) and I suspect that up to 50 at most - try magnifying the plot in the interesting region to see the details. This implies that the sampling of the data is just too dense - you don't get any new information but you make huge additional effort to achieve, in fact, frequencies with power practically equal to zero. The approach in image 3 - cutting the data in half - does not change anything because you keep the same sampling step. In images 4 & 5 you: (i) made the sampling significantly more sparse, and (ii) added noise, i.e. added several frequencies not present in the original data set. That's why those spectra look "nicer" - they just have more details to show. Finally, due to finiteness of the data set you'll get ambiguous recognitions of frequencies due to folding those that are outside the Nyquist interval into this interval. That's why I suggest to use the Lomb-Scargle periodogram - it tests only the pre-defined periods. Moreover, you could also use a wavelet scalogram to investigate which frequencies are present at which times in the time series. To conclude, your code seems to be ok (at least at first glance); I'd suggest digging deeper into the setting and mathematical features of the FFT algorithm to grasp its limitations. But that is not a Mathematica related question anymore.


You might want to use a Lomb-Scargle periodogram: generate some data:

t = Table[i, {i, 30., 60., 0.025}];
x = Table[Sin[i] Sin[2 i], {i, 30., 60., 0.025}];
tx = Table[{t[[i]], x[[i]]}, {i, 1, Length[t]}];
xbar = Mean[x];
std = 1/(Length[x] - 1.) Sum[(x[[i]] - xbar)^2, {i, 1, Length[x]}];

enter image description here

Definition of the LS periodogram:

LS[tx_, nP_] := 
 LS[tx, nP] = 
  Block[{t, x, xbar, std, \[Tau], P, \[CapitalOmega], diff, lim, 
     periodogram},
    t = Table[tx[[i, 1]], {i, 1, Length[tx]}];
    x = Table[tx[[i, 2]], {i, 1, Length[tx]}];
    xbar = Mean[x];
    std = 1/(Length[x] - 1.) Sum[(x[[i]] - xbar)^2, {i, 1, Length[x]}];
    \[Tau][\[Omega]_] := \[Tau][\[Omega]] = ArcTan[\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(Length[tx]\)]\(Sin[
         2. \ \[Omega]\ x[[i]]]\)\)/\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(Length[tx]\)]\(Cos[
         2. \ \[Omega]\ x[[i]]]\)\)]/(2. \[Omega]);
    P[\[Omega]_] := P[\[Omega]] = (0.5/std) ((\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(Length[
             tx]\)]\(\((x[[i]] - 
              xbar)\)\ Cos[\[Omega]\ \((t[[
                 i]] - \[Tau][\[Omega]])\)]\)\))^2/\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(Length[tx]\)]
\*SuperscriptBox[\(Cos[\[Omega]\ \((t[[
                i]] - \[Tau][\[Omega]])\)]\), \(2\)]\) + (\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(Length[

             tx]\)]\(\((x[[i]] - 
              xbar)\)\ Sin[\[Omega]\ \((t[[
                 i]] - \[Tau][\[Omega]])\)]\)\))^2/\!\(
\*UnderoverscriptBox[\(\[Sum]\), \(i = 1\), \(Length[tx]\)]
\*SuperscriptBox[\(Sin[\[Omega]\ \((t[[
                i]] - \[Tau][\[Omega]])\)]\), \(2\)]\));
    \[CapitalOmega][T_] := \[CapitalOmega][T] = 2. \[Pi]/T;
    diff = Table[t[[i + 1]] - t[[i]], {i, 1, Length[t] - 1}]; 
    lim = \[CapitalOmega] /@ {Last[t] - First[t], 2 Min[diff]};
    periodogram = 
     ParallelTable[{\[Omega], P[\[Omega]]}, {\[Omega], lim[[1]], 
       lim[[2]], (lim[[2]] - lim[[1]])/nP}]
    ]

The argumens of the function are tx_ = input data (time series) of the form $(t, x)$ and nP_ = number of frequencies/periods to be tested; the function returns par as $(frequency, its\, power)$. Let's set

nP=12010;

Then

perdogram = LS[tx, nP]; // AbsoluteTiming

took about 20 seconds. Finally

ListPlot[{Take[perdogram, 400], Take[perdogram, 400]}, Frame -> True, 
 PlotStyle -> {Black, Red}, Joined -> {True, False}, PlotRange -> All,
  ImageSize -> 500, 
 FrameLabel -> {"\[Omega]", Rotate["P(\[Omega])", 270 Degree]}, 
 AxesOrigin -> {0, 0}]

enter image description here

To get $p$-values for a specific frequency, define:

T[periodogram_, rankmax_] := 
 T[periodogram, rankmax] = Block[{power, \[Omega]max, period},
   power = Table[periodogram[[i, 2]], {i, 1, Length[periodogram]}];
   \[Omega]max = 
    periodogram[[Position[power, RankedMax[power, rankmax]][[1, 1]], 
     1]];
   {period = 2 \[Pi]/\[Omega]max, RankedMax[power, rankmax]}
   ]

where its arguments are: the periodogram_, i.e. the output of LS[tx_, nP_], and rankmax_, i.e. the number of a maximum in the distribution of power of a periodogram; it returns the period and its power.

For example:

T[perdogram, 1]
T[perdogram, 2]

give

{6.32827, 301.723}
{2.10279, 294.189}

The last step is to define a $p$-value as

p[z_, n_] := p[z, n] = 1 - (1 - Exp[-z])^n

where z_ is the value of the power at a given frequency/period (i.e., at a peak), and n_ is the number of frequencies/periods testes (i.e., it's nP_ from te definition of LS). E.g.,

p[301.723, nP]
p[294.189, nP]
p[16.066488, nP]
p[12.652132, nP]

gives

0.
0.
0.00126381
0.0377113

You can also plot the function for the $p$-value:

Plot[p[z, nP], {z, 0, 20}, Frame -> True, PlotStyle -> Black, FrameLabel -> {"z", "p-value"}]

enter image description here

$\endgroup$
  • $\begingroup$ thanks for the answer, still I don't get why I have a vertical line in the origin and zero everywhere else, like in image 2 and 3. $\endgroup$ – mtlvc0 Aug 15 '16 at 20:08
  • $\begingroup$ As stated above in the comments, in FFT you are mainly interested in the range $(0,5000)$ (referring to your image 2); judging from the graph in image 1 (which is highly periodic - the damping doesn't matter much) one should expect a few, relatively low, frequencies: the lowest corresponding to the length of the time series (hence not informative) and the others related to the frequencies indeed present in the data set (magnifiy the plot to reveal te details). In general, you'll have a problem with identifying low- and high-frequency terms, that's why I suggest to use an LS periodogram. $\endgroup$ – corey979 Aug 15 '16 at 23:05
  • $\begingroup$ The notes from the comment of @Hugh above are also a good reference. $\endgroup$ – corey979 Aug 15 '16 at 23:26
  • $\begingroup$ thanks a lot, i checked those notes but there weren't some concepts that I needed, I will magnify and see and I will use the periodogram, it's the first time I hear about that and I couldn't understand the intrinsic advantage $\endgroup$ – mtlvc0 Aug 16 '16 at 8:20

Not the answer you're looking for? Browse other questions tagged or ask your own question.