intfft

INTFFT - calculate the tune from interpolated FFT of the trajectory.
 INTFFT(X) X must be a column vector.
  If X is a matrix - each column is treated as
  a separate trajectory
 INTFFT(X,GUESS,DELTA) searches for peaks in the FFT spectrum
  only within the range (X-DELTA ... X+DELTA) 
  The same values of GUESS and DELTA are used for all columns of X

0001 function tune = intfft(X,varargin);
0002 %INTFFT - calculate the tune from interpolated FFT of the trajectory.
0003 % INTFFT(X) X must be a column vector.
0004 %  If X is a matrix - each column is treated as
0005 %  a separate trajectory
0006 % INTFFT(X,GUESS,DELTA) searches for peaks in the FFT spectrum
0007 %  only within the range (X-DELTA ... X+DELTA)
0008 %  The same values of GUESS and DELTA are used for all columns of X
0009 
0010 
0011 [N,L] = size(X);
0012 if L == 0;
0013     tune = NaN;
0014     
0015     return
0016 end
0017 % apply hanning window
0018 %W = diag(sin(pi*(0:N-1)/(N-1)).^2);
0019 
0020 %XFFTABS = abs(fft(W*X));
0021 XFFTABS = abs(fft(X));
0022 %Z = zeros(size(XFFTABS));
0023 if nargin==3
0024     GUESS = varargin{1};
0025     DELTA = varargin{2};
0026 %     LR = floor(N*(GUESS-DELTA));
0027 %     UR = ceil(N*(GUESS+DELTA));
0028 %     Z(sub2ind(size(XFFTABS),LR,1:length(LR))) = 1;
0029 %     Z(sub2ind(size(XFFTABS),UR+1,1:length(UR)))= -1;
0030     
0031     
0032     searchrange = floor(N*(GUESS-DELTA)):ceil(N*(GUESS+DELTA));
0033     
0034     %[psi_k,k] = max(cumsum(Z).*XFFTABS);
0035     [psi_k,k] = max(XFFTABS(searchrange,:));
0036     k=k+floor(N*(GUESS-DELTA))-1; 
0037     
0038 else 
0039     [psi_k,k] = max(XFFTABS(1:floor(N/2),:));
0040 end
0041 
0042 psi_k_plus  = XFFTABS((k+1)+(0:(L-1))*N);
0043 psi_k_minus = XFFTABS((k-1)+(0:(L-1))*N);
0044 
0045 G = psi_k_plus>psi_k_minus;
0046 
0047 k_r = k+G;
0048 k_l = k-~G;
0049 
0050 psi_l = XFFTABS(k_l+(0:(L-1))*N);
0051 psi_r = XFFTABS(k_r+(0:(L-1))*N);
0052 
0053 tune = (k_l-1)/N + atan( psi_r*sin(pi/N)./(psi_l + psi_r*cos(pi/N)))/pi;