home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / POLYFUN.DI$ / INTERPFT.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.2 KB  |  51 lines
  1. function y = interpft(x,ny);
  2. %INTERPFT 1-D interpolation using a FFT method.
  3. %    Y = INTERPFT(X,N) returns a length N vector Y obtained by
  4. %    interpolation in the Fourier transform of X.  N must be larger
  5. %    than the length of X.  The result has circular periodicity.
  6. %
  7. %    See also INTERP1, INTERP2, GRIDDATA.
  8.  
  9. %    Charles R. Denham, MathWorks, 1988.
  10. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  11.  
  12. if nargin < 1      % Demonstration.
  13.    nx = 10;
  14.    x= rand(nx,1);
  15.    ny = 100;
  16. end
  17.  
  18. [m,n] = size(x);
  19. oldm = m;
  20. if m == 1
  21.    x = x.';
  22.    [m,n] = size(x);
  23. end
  24. if ny <= m, error('NY must exceed the length of X.'), end
  25. i = m/2+1;
  26. if rem(m,2) ~= 0, i = (m+1)/2; end
  27. a = fft(x);
  28. b = [a(1:i,:) ; zeros(ny-m,n) ; a(i+1:m,:)];
  29. if rem(m,2) == 0
  30.    if rem(ny,2) == 0
  31.       b(i) = b(i)/2;
  32.       b(i+ny-m) = b(i);
  33.      else
  34.       b(i)=0;
  35.    end
  36. end
  37. y = ifft(b);
  38. if oldm == 1, y = y.'; end
  39. y = y * ny / m;
  40. if ~any(imag(x)), y = real(y); end
  41.  
  42. if nargin < 1      % Demonstration.
  43.    hold off
  44.    ax = (0:(ny-1)) ./ ny; ax = nx .* ax + 1;
  45.    plot(ax, y), hold on
  46.    ax = (0:(nx-1)) ./ nx; ax = nx .* ax + 1;
  47.    plot(ax, x, '*'), hold off
  48.    title('Demonstration of INTERPFT')
  49.    xlabel('x'), ylabel('y')
  50. end
  51.