octave_packages/m/general/interpft.m

   1 ## Copyright (C) 2001-2012 Paul Kienzle
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} interpft (@var{x}, @var{n})
  21 ## @deftypefnx {Function File} {} interpft (@var{x}, @var{n}, @var{dim})
  22 ##
  23 ## Fourier interpolation.  If @var{x} is a vector, then @var{x} is
  24 ## resampled with @var{n} points.  The data in @var{x} is assumed to be
  25 ## equispaced.  If @var{x} is an array, then operate along each column of
  26 ## the array separately.  If @var{dim} is specified, then interpolate
  27 ## along the dimension @var{dim}.
  28 ##
  29 ## @code{interpft} assumes that the interpolated function is periodic,
  30 ## and so assumptions are made about the endpoints of the interpolation.
  31 ##
  32 ## @seealso{interp1}
  33 ## @end deftypefn
  34
  35 ## Author: Paul Kienzle
  36 ## 2001-02-11
  37 ##    * initial version
  38 ## 2002-03-17 aadler
  39 ##    * added code to work on matrices as well
  40 ## 2006-05-25 dbateman
  41 ##    * Make it matlab compatiable, cutting out the 2-D interpolation
  42
  43 function z = interpft (x, n, dim)
  44
  45   if (nargin < 2 || nargin > 3)
  46     print_usage ();
  47   endif
  48
  49   if (! (isscalar (n) && n == fix (n)))
  50     error ("interpft: N must be a scalar integer");
  51   endif
  52
  53   if (nargin == 2)
  54     if (isrow (x))
  55       dim = 2;
  56     else
  57       dim = 1;
  58     endif
  59   endif
  60
  61   nd = ndims (x);
  62
  63   if (dim < 1 || dim > nd)
  64     error ("interpft: invalid dimension DIM");
  65   endif
  66
  67   perm = [dim:nd, 1:(dim-1)];
  68   x = permute (x, perm);
  69   m = rows (x);
  70
  71   inc = max (1, fix (m/n));
  72   y = fft (x) / m;
  73   k = floor (m / 2);
  74   sz = size (x);
  75   sz(1) = n * inc - m;
  76
  77   idx = repmat ({':'}, nd, 1);
  78   idx{1} = 1:k;
  79   z = cat (1, y(idx{:}), zeros (sz));
  80   idx{1} = k+1:m;
  81   z = cat (1, z, y(idx{:}));
  82   z = n * ifft (z);
  83
  84   if (inc != 1)
  85     sz(1) = n;
  86     z = inc * reshape (z(1:inc:end), sz);
  87   endif
  88
  89   z = ipermute (z, perm);
  90
  91 endfunction
  92
  93
  94 %!demo
  95 %! t = 0 : 0.3 : pi; dt = t(2)-t(1);
  96 %! n = length (t); k = 100;
  97 %! ti = t(1) + [0 : k-1]*dt*n/k;
  98 %! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
  99 %! yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
 100 %! plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ...
 101 %!       ti, interpft (y, k), 'c', t, y, 'r+');
 102 %! legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data');
 103
 104 %!shared n,y
 105 %! x = [0:10]'; y = sin(x); n = length (x);
 106 %!assert (interpft(y, n), y, 20*eps);
 107 %!assert (interpft(y', n), y', 20*eps);
 108 %!assert (interpft([y,y],n), [y,y], 20*eps);
 109
 110 %% Test input validation
 111 %!error interpft ()
 112 %!error interpft (1)
 113 %!error interpft (1,2,3)
 114 %!error (interpft(1,[n,n]))
 115 %!error (interpft(1,2,0))
 116 %!error (interpft(1,2,3))