Message Boards Message Boards

1
|
5658 Views
|
5 Replies
|
4 Total Likes
View groups...
Share
Share this post:

Signal or (quasi periodic) trajectory frequency in Mathematica

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 convolution
on 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.
POSTED BY: Ivan Morozov
5 Replies
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 BY: Ivan Morozov
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]
POSTED BY: Ivan Morozov
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_13

I.M.
POSTED BY: Ivan Morozov
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 BY: Ivan Morozov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract