calcnaff

 [nu amp phase] = calcnaff(Y, Yp, Win)

  INPUTS
  1. Y  - position vector
  2. Yp -
  3. WindowType  - Window type - 0 {Default} no windowing
                                 1 Window of Hann
                                 2 etc
  4. nfreq - Maximum number of fundamental frequencies to search for
             10 {Default}
  5. debug - 1 means debug flag turned on
             0 {Default}

  Optional Flags
  'Debug' - turn on deubbing flag
  'Display' - print ou results
  'Hanning' - select Window of Hann, WindowType = 1
  'Raw' or 'NoWindow' - select Window of Hann, WindowType = 0

  OUTPUTS
  1. frequency - frequency vector with sorted amplitudes
                 by default the algorithm try to compute the 10 first fundamental
                 frequencies of the system.
  2. amplitude - amplitudes associated with fundamental frequencies
  3. phase - phases associated with fundamental frequencies

  NOTES
  1. Mimimum number of turns is 64 (66)
  2. Number of turn has to be a multiple of 6 for internal optimization
  reason and just above a power of 2. Example 1026 is divived by 6 and
  above 1024 = pow2(10)

  Examples
  NT = 9996; % divided by 6
  simple quasiperiodic (even period) motion 
  y =2+0.1*cos(pi*(0:NT-1))+0.00125*cos(pi/3*(0:NT-1));
  yp=2+0.1*sin(pi*(0:NT-1))+0.00125*sin(pi/3*(0:NT-1));
 
  [nu ampl phase] = calcnaff(y,yp,1); % with windowing

0001 function [frequency,amplitude,phase] = calcnaff(Y, Yp, varargin)
0002 % [nu amp phase] = calcnaff(Y, Yp, Win)
0003 %
0004 %  INPUTS
0005 %  1. Y  - position vector
0006 %  2. Yp -
0007 %  3. WindowType  - Window type - 0 {Default} no windowing
0008 %                                 1 Window of Hann
0009 %                                 2 etc
0010 %  4. nfreq - Maximum number of fundamental frequencies to search for
0011 %             10 {Default}
0012 %  5. debug - 1 means debug flag turned on
0013 %             0 {Default}
0014 %
0015 %  Optional Flags
0016 %  'Debug' - turn on deubbing flag
0017 %  'Display' - print ou results
0018 %  'Hanning' - select Window of Hann, WindowType = 1
0019 %  'Raw' or 'NoWindow' - select Window of Hann, WindowType = 0
0020 %
0021 %  OUTPUTS
0022 %  1. frequency - frequency vector with sorted amplitudes
0023 %                 by default the algorithm try to compute the 10 first fundamental
0024 %                 frequencies of the system.
0025 %  2. amplitude - amplitudes associated with fundamental frequencies
0026 %  3. phase - phases associated with fundamental frequencies
0027 %
0028 %  NOTES
0029 %  1. Mimimum number of turns is 64 (66)
0030 %  2. Number of turn has to be a multiple of 6 for internal optimization
0031 %  reason and just above a power of 2. Example 1026 is divived by 6 and
0032 %  above 1024 = pow2(10)
0033 %
0034 %  Examples
0035 %  NT = 9996; % divided by 6
0036 %  simple quasiperiodic (even period) motion
0037 %  y =2+0.1*cos(pi*(0:NT-1))+0.00125*cos(pi/3*(0:NT-1));
0038 %  yp=2+0.1*sin(pi*(0:NT-1))+0.00125*sin(pi/3*(0:NT-1));
0039 %
0040 %  [nu ampl phase] = calcnaff(y,yp,1); % with windowing
0041 
0042 
0043 % Written by Laurent S. Nadolski
0044 % April 6th, 2007
0045 % Modification September 2009:
0046 %  test if constant data or nan data
0047 
0048 % BUG in nafflib: returns nan even if valid data. Number of try
0049 nitermax = 10;
0050 
0051 % Flag factory
0052 [wraw1,args]=getflag(varargin,'Raw'); %#ok<ASGLU>
0053 [wraw2,args]=getflag(args,'NoWindow'); %#ok<ASGLU>
0054 [whann,args]=getflag(args,'Hanning');
0055 [dbg,args]=getflag(args,'Debug');
0056 [DisplayFlag,args]=getflag(args,'Display');
0057 [WindowType,nfreq,DebugFlag]=getargs(args,{0,10,double(dbg)});
0058 if whann, WindowType=1; end
0059 
0060 
0061 % Test wether nan or constant data
0062 if any(isnan(Y(1,:)))
0063     fprintf('Warning Y contains NaNs\n');
0064     frequency =NaN; amplitude = NaN;  phase = NaN;
0065 elseif any(isnan(Y(1,:)))
0066     fprintf('Warning Yp contains NaNs\n');
0067     frequency =NaN; amplitude = NaN;  phase = NaN;
0068 elseif (mean(Y) == Y(1) && mean(Yp) == Yp(1))
0069     fprintf('Warning data are constant\n');
0070     frequency = 0; amplitude = 0;  phase = 0;
0071 else % Frequency map analysis
0072     [frequency,amplitude,phase] = nafflib(Y, Yp, WindowType,nfreq,DebugFlag);
0073     %It seems there is a bug in nafflib, something returns nan even for valid data
0074     niter = 0;
0075     while any(isnan(frequency)) && (niter < nitermax)
0076         pause(2);
0077         fprintf('Warning Nan returned by NAFF (x%d)\n', niter);
0078         niter = niter +1;
0079         [frequency,amplitude,phase] = nafflib(Y, Yp, WindowType,nfreq,1); % add debugging
0080     end
0081         
0082     if DisplayFlag
0083         fprintf('*** Naff run on %s\n', datestr(clock))
0084         for k = 1:length(frequency)
0085             fprintf('nu(%d) =% .15f amplitude = % .15f phase = % .15f \n', ...
0086                 k,frequency(k), amplitude(k),phase(k));
0087         end
0088     end
0089 end