function retval = FFTsort(p)
% usage: ndx = FFTsort(p)
% description: given positive integer p as input,
% returns an index vector of length N = 2^p
% suitable for implementing the FFT, with
% indexing starting at 1.
 
% programmer: t. shores
% last update: 2/6/16
 
% Octave and Matlab compatible file
N = 2^p;
retval = (1:N);
stride = N;
for k = p:-1:1
  for j = 1:stride:N-1
    m = j+stride-1;
    retval(j:m) = [retval(j:2:m),retval(j+1:2:m)];
  end
  stride = stride/2;
end
