function y = interpft(x,ny);
%INTERPFT 1-D interpolation using a FFT method.
%	Y = INTERPFT(X,N) returns a length N vector Y obtained by
%	interpolation in the Fourier transform of X.  N must be larger
%	than the length of X.  The result has circular periodicity.
%
%	See also INTERP1, INTERP2, GRIDDATA.

%	Charles R. Denham, MathWorks, 1988.
%	Copyright (c) 1984-93 by The MathWorks, Inc.

if nargin < 1      % Demonstration.
   nx = 10;
   x= rand(nx,1);
   ny = 100;
end

[m,n] = size(x);
oldm = m;
if m == 1
   x = x.';
   [m,n] = size(x);
end
if ny <= m, error('NY must exceed the length of X.'), end
i = m/2+1;
if rem(m,2) ~= 0, i = (m+1)/2; end
a = fft(x);
b = [a(1:i,:) ; zeros(ny-m,n) ; a(i+1:m,:)];
if rem(m,2) == 0
   if rem(ny,2) == 0
      b(i) = b(i)/2;
      b(i+ny-m) = b(i);
     else
      b(i)=0;
   end
end
y = ifft(b);
if oldm == 1, y = y.'; end
y = y * ny / m;
if ~any(imag(x)), y = real(y); end

if nargin < 1      % Demonstration.
   hold off
   ax = (0:(ny-1)) ./ ny; ax = nx .* ax + 1;
   plot(ax, y), hold on
   ax = (0:(nx-1)) ./ nx; ax = nx .* ax + 1;
   plot(ax, x, '*'), hold off
   title('Demonstration of INTERPFT')
   xlabel('x'), ylabel('y')
end
