# Signal or (quasi periodic) trajectory frequency in Mathematica

Posted 10 years ago
5172 Views
|
5 Replies
|
4 Total Likes
|
 Hi,Continuing the topic started in here I've decided to share my frequency search routine.Note that I do not have Newton algorithm in it and the routing is not optimized.The idea behind adding Newton algorithm is the following: when first approximation of the frequency is found instead of using convolutionon the full signal convolution is performed for "left" and "right" regions separately. The best is then choosen and division is performed again.Division should be performed until some pre-defined accuracy for frequency is reached (typically ~10^-9...10^-10)This module can be used only for equally spaced data (in time or space).If someone is interested in improving this code I'll be happy to assist. Unfortunately, I have no time to finish it by myself.I.M.
5 Replies
Sort By:
Posted 10 years ago
 Hi,here is an updated version for frequency identification. I've cleaned up the code and added iterative search.Now the accuracy can be near machine precision with relatively low sample size. NAFF[ARG_List] := Module[    {     (* LOCALS *)     LENGTH,     WINDOW,     FOURIER,     POINT1, POINT2, POINT3,     OMEGA1, OMEGA2, OMEGA3,     ORBIT    },   LENGTH = Length@ARG;   WINDOW = Table[0.5 (1. - Cos[2 Pi j/(LENGTH - 1)]), {j, 1, LENGTH}];   ORBIT = (ARG - Mean[ARG]) WINDOW;   FOURIER = Abs@Fourier[ORBIT];   POINT1 = Position[FOURIER, Max[FOURIER]][[1, 1]];   OMEGA1 = N[(POINT1 - 2)/LENGTH];   ORBIT = ORBIT  Exp[2 Pi I OMEGA1 N[Range[0, LENGTH - 1]]] WINDOW;   FOURIER = Abs@Fourier[ORBIT, FourierParameters -> {0, 2/LENGTH}];   POINT2 = Position[FOURIER, Max[FOURIER]][[1, 1]];   OMEGA2 = N[1/LENGTH (POINT1 - 2 + 2 (POINT2 - 1)/LENGTH)];   Off[FindMaximum::lstol];   POINT3 = FindMaximum[Interpolation[FOURIER, InterpolationOrder -> 8][s], {s,POINT2}, Method -> "Newton",WorkingPrecision -> 20][[2, 1, 2]];   OMEGA3 = N[1/LENGTH (POINT1 - 2 + 2 (POINT3 - 1)/LENGTH)];   OMEGA3   ];NAFF::usage = "Takes quasi-periodic signal as input  and returnes (1st fundumental) frequncy (i.e. with largest amplitude)" ;(* EG-1 *)f1[x_] := 0.8 Sin[21 Pi/E x] + 0.75 Sin[13 Pi/E x];data = Table[f1[x], {x, 0, 2048, 1}];v = 1. - N@Mod[21 Pi/E, 2 Pi]/(2 Pi)NAFF[data]Abs[%% - %](* EG-2 *)f2[x_] := 0.8 Sin[21 Pi/E x] + 0.85 Sin[13 Pi/E x];data = Table[f2[x], {x, 0, 2048, 1}];v = N@Mod[13 Pi/E, 2 Pi]/(2 Pi)NAFF[data]Abs[%% - %]I.M.
Posted 10 years ago
 Hi, Daniel,>> Maybe similar to what shows up in this MSE post?Quite similar indeed, but at that post the emphasis is on speed, here -- on accuracy.Some applications require as much accuracy as one can get, for example,stable vs. unstable (chaotic) motion for Hamiltonian systems, e.g.:http://link.springer.com/chapter/10.1007%2F978-94-011-4673-9_13I.M.
Posted 10 years ago
 Maybe similar to what shows up in this MSE post?http://mathematica.stackexchange.com/questions/42783/how-do-you-find-the-frequency-and-amplitude-from-fourier
Posted 10 years ago
 CODE:  FrequencySearch[data_, step_, window_, test_: {"skip"}] := Module[     {      DATA, (* EQUALY SPACED SIGNAL *)      WINDOW, (* WINDOW FUNCTION *)              POINTS, (* # OF POINTS IN DATA *)      STEP, (*       POINTS SPACING *)     WINDOWFUNC, X, WINDOWDATA, FOURIER, (*      LOCALS *)     POS1, POS2, POS3,(* LOCALS *)     T1, T2, (*      PERIOD WITHOUT SIGNAL CONVOLUTION AND WITH IT *)     OMEGA1,      OMEGA2, (* SAME STAFF FOR FREQUENCY *)     TEST, (* A DEBUG FLAG *)           POSREQ, OMEGAREQ, TREQ, RUN (* LOCALS *)     },    {     (* ASSIGN INPUT TO LOCAL VARS *)     TEST = test;     STEP = step;     DATA = data;     POINTS = Length[DATA];     WINDOW = window;                  (* TESTING*)             If[TEST == "Test", {                   Print["STEP 1: PLOT INITIAL DATA"];                   Print[ListPlot[DATA]];               }];          (* REMOVE SIGNAL AVERAGE *)     DATA = DATA - Mean[DATA];                          (* TESTING*)             If[TEST == "Test", {                   Print["STEP 2: REMOVE SIGNAL AVERAGE VALUE"];                   Print[ListPlot[DATA]];               }];          (* APPLY WINDOW *)          If[WINDOW ==        "Hann", {WINDOWFUNC[X_] :=          0.5 (1 - Cos[2 Pi X/(POINTS - 1)]);}, {WINDOWFUNC[X_] :=          0.42 - 0.5 Cos[2 Pi X/(POINTS - 1)] +           0.08 Cos[4 Pi X/(POINTS - 1)];}];     WINDOWDATA = Table[WINDOWFUNC[X], {X, 1, POINTS}];     DATA = DATA WINDOWDATA;                          (* TESTING*)             If[TEST == "Test", {                   Print["STEP 3: APPLY WINDOW FUNCTION"];                   Print[ListPlot[DATA]];               }];          (* FOURIER WITHOUT CONVOLUTION *)     FOURIER = Abs[Fourier[DATA]];     POS1 = Position[FOURIER, Max[FOURIER]][[1, 1]];     T1 = POINTS/(POS1 - 2 ) STEP;     OMEGA1 = (POS1 - 2)/POINTS;                          (* TESTING*)             If[TEST == "Test", {       Print["STEP 4: FOURIER TRANSFORM AND LOCATE INITIAL MAXIMUM"];                   Print[Show[                         ListPlot[FOURIER],                         Graphics[{Red, Point[{POS1, FOURIER[[POS1]]}]}],          PlotRange -> All                     ]];                   Print["Initial position: ", POS1];                   Print["Initial frequency: ", 2 Pi/T1 // N];               }];          (* OVERLAP/CONVOLUTION *)          DATA = DATA Exp[2 Pi I OMEGA1 N[Range[0, POINTS - 1]]];                          (* TESTING*)             If[TEST == "Test", {                   Print["STEP 5: OVERLAP (Re,Im,Abs)"];                   Print[ListPlot[Re[DATA]]];                   Print[ListPlot[Im[DATA]]];                   Print[ListPlot[Abs[DATA]]];               }];          (* FOURIER WITH CONVOLUTION *)          FOURIER = Abs[Fourier[DATA, FourierParameters -> {0, 2/POINTS}]];     POS2 = Position[FOURIER, Max[FOURIER]][[1, 1]];     T2 = POINTS/(POS1 - 2 + 2 (POS2 - 1)/POINTS) STEP;             (* TESTING*)             If[TEST == "Test", {               Print["STEP 6: FOURIER TRANSFORM AND LOCATE NEW MAXIMUM"];                   Print[Show[                        ListPlot[FOURIER],                        Graphics[{Red, Point[{POS2, FOURIER[[POS2]]}]}],         PlotRange -> All                    ]];                  Print["New position: ", POS2];                  Print["New frequency: ", 2 Pi/T2 // N];              }];    (* RESULTS *)    OMEGA1 = 2 Pi/T1 (* NO CONVOLUTION *);    OMEGA2 = 2 Pi /T2(* CONVOLUTION *);    {N[OMEGA1], N[OMEGA2]}    }];
Posted 10 years ago
 EXAMPLE (*TEST*) f[x_] := 1 + Sin[21 Pi/E x]; fr = N[21 Pi/E]; data = Table[f[x], {x, 0, 5, 0.0001}]; FrequencySearch[data, 0.0001, "Hann", "Test"]; (* NO CONVOLUTION *) f1 = %[[1, 1]] (* CONVOLUTION *) f2 = %%[[1, 2]](* EXACT *)fr(* ERROR *)Abs[(f1 - fr)/fr]Abs[(f2 - fr)/fr]
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments