1 | function [frequency amplitude phase] = calcnaff(Y, Yp, varargin) |
---|
2 | % [nu amp phase] = calcnaff(Y, Yp, Win) |
---|
3 | % |
---|
4 | % INPUTS |
---|
5 | % 1. Y - position vector |
---|
6 | % 2. Yp - |
---|
7 | % 3. WindowType - Window type - 0 {Default} no windowing |
---|
8 | % 1 Window of Hann |
---|
9 | % 2 etc |
---|
10 | % 4. nfreq - Maximum number of fundamental frequencies to search for |
---|
11 | % 10 {Default} |
---|
12 | % 5. debug - 1 means debug flag turned on |
---|
13 | % 0 {Default} |
---|
14 | % |
---|
15 | % Optional Flags |
---|
16 | % 'Debug' - turn on deubbing flag |
---|
17 | % 'Display' - print ou results |
---|
18 | % 'Hanning' - select Window of Hann, WindowType = 1 |
---|
19 | % 'Raw' or 'NoWindow' - select Window of Hann, WindowType = 0 |
---|
20 | % |
---|
21 | % OUTPUTS |
---|
22 | % 1. frequency - frequency vector with sorted amplitudes |
---|
23 | % by default the algorithm try to compute the 10 first fundamental |
---|
24 | % frequencies of the system. |
---|
25 | % 2. amplitude - amplitudes associated with fundamental frequencies |
---|
26 | % 3. phase - phases associated with fundamental frequencies |
---|
27 | % |
---|
28 | % NOTES |
---|
29 | % 1. Mimimum number of turns is 64 (66) |
---|
30 | % 2. Number of turn has to be a multiple of 6 for internal optimization |
---|
31 | % reason and just above a power of 2. Example 1026 is divived by 6 and |
---|
32 | % above 1024 = pow2(10) |
---|
33 | % |
---|
34 | % Examples |
---|
35 | % NT = 9996; % divided by 6 |
---|
36 | % simple quasiperiodic (even period) motion |
---|
37 | % y =2+0.1*cos(pi*(0:NT-1))+0.00125*cos(pi/3*(0:NT-1)); |
---|
38 | % yp=2+0.1*sin(pi*(0:NT-1))+0.00125*sin(pi/3*(0:NT-1)); |
---|
39 | % |
---|
40 | % [nu ampl phase] = calcnaff(y,yp,1); % with windowing |
---|
41 | |
---|
42 | |
---|
43 | % Written by Laurent S. Nadolski |
---|
44 | % April 6th, 2007 |
---|
45 | % Modification September 2009: |
---|
46 | % test if constant data or nan data |
---|
47 | |
---|
48 | % Default flags |
---|
49 | DebugFlag = 0; % 0/1 dubugging flag in NAFF |
---|
50 | WindowType = 1; % Window type |
---|
51 | nfreq = 10; % number of frequency |
---|
52 | DisplayFlag = 0; % Display Figure |
---|
53 | |
---|
54 | % BUG in nafflib: returns nan even if valid data. Number of try |
---|
55 | nitermax = 10; |
---|
56 | |
---|
57 | % Flag factory |
---|
58 | for ik = length(varargin):-1:1 |
---|
59 | if strcmpi(varargin{ik},'Debug') |
---|
60 | DebugFlag = 1; |
---|
61 | varargin(ik) = []; |
---|
62 | elseif strcmpi(varargin{ik},'NoDebug') |
---|
63 | DebugFlag = 0; |
---|
64 | varargin(ik) = []; |
---|
65 | elseif strcmpi(varargin{ik},'Display') |
---|
66 | DisplayFlag = 1; |
---|
67 | varargin(ik) = []; |
---|
68 | elseif strcmpi(varargin{ik},'NoDisplay') |
---|
69 | DisplayFlag = 0; |
---|
70 | varargin(ik) = []; |
---|
71 | elseif strcmpi(varargin{ik},'Hanning') |
---|
72 | WindowType = 1; |
---|
73 | varargin(ik) = []; |
---|
74 | elseif strcmpi(varargin{ik},'NoWindow') || strcmpi(varargin{ik},'Raw') |
---|
75 | WindowType = 1; |
---|
76 | varargin(ik) = []; |
---|
77 | end |
---|
78 | end |
---|
79 | |
---|
80 | if length(varargin) >= 1 |
---|
81 | WindowType = varargin{1}; |
---|
82 | end |
---|
83 | |
---|
84 | if length(varargin) >= 2 |
---|
85 | nfreq = varargin{2}; |
---|
86 | end |
---|
87 | |
---|
88 | if length(varargin) >= 3 |
---|
89 | DebugFlag = varargin{3}; |
---|
90 | end |
---|
91 | |
---|
92 | if length(varargin) >= 4 |
---|
93 | error('Too many arguments'); |
---|
94 | end |
---|
95 | |
---|
96 | % Test wether nan or constant data |
---|
97 | if any(isnan(Y(1,:))) |
---|
98 | fprintf('Warning Y contains NaNs\n'); |
---|
99 | frequency =NaN; amplitude = NaN; phase = NaN; |
---|
100 | elseif any(isnan(Y(1,:))) |
---|
101 | fprintf('Warning Yp contains NaNs\n'); |
---|
102 | frequency =NaN; amplitude = NaN; phase = NaN; |
---|
103 | elseif (mean(Y) == Y(1) && mean(Yp) == Yp(1)) |
---|
104 | fprintf('Warning data are constant\n'); |
---|
105 | frequency = 0; amplitude = 0; phase = 0; |
---|
106 | else % Frequency map analysis |
---|
107 | [frequency amplitude phase] = nafflib(Y, Yp, WindowType,nfreq,DebugFlag); |
---|
108 | %It seems there is a bug in nafflib, something returns nan even for valid data |
---|
109 | niter = 0; |
---|
110 | while any(isnan(frequency)) && (niter < nitermax) |
---|
111 | pause(2); |
---|
112 | fprintf('Warning Nan returned by NAFF (x%d)\n', niter); |
---|
113 | niter = niter +1; |
---|
114 | [frequency amplitude phase] = nafflib(Y, Yp, WindowType,nfreq,1); % add debugging |
---|
115 | end |
---|
116 | |
---|
117 | if DisplayFlag |
---|
118 | fprintf('*** Naff run on %s\n', datestr(clock)) |
---|
119 | for k = 1:length(frequency) |
---|
120 | fprintf('nu(%d) =% .15f amplitude = % .15f phase = % .15f \n', ... |
---|
121 | k,frequency(k), amplitude(k),phase(k)); |
---|
122 | end |
---|
123 | end |
---|
124 | end |
---|
125 | |
---|