← Back to branch summary

~ubuntu-branches/debian/sid/octave3.0/sid

« back to all changes in this revision

Viewing changes to scripts/general/interpft.m

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2007-12-23 16:04:15 UTC
  • Revision ID: james.westby@ubuntu.com-20071223160415-n4gk468dihy22e9v
Tags: upstream-3.0.0
Import upstream version 3.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/general/interpft.m
 
1
## Copyright (C) 2001, 2006, 2007 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 end points 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 (nargin == 2)
 
50
    if (isvector (x) && size (x, 1) == 1)
 
51
      dim = 2;
 
52
    else
 
53
      dim = 1;
 
54
    endif
 
55
  endif
 
56
 
 
57
  if (! isscalar (n))
 
58
    error ("interpft: n must be an integer scalar");
 
59
  endif
 
60
 
 
61
  nd = ndims (x);
 
62
 
 
63
  if (dim < 1 || dim > nd)
 
64
    error ("interpft: integrating over invalid dimension");
 
65
  endif
 
66
 
 
67
  perm = [dim:nd, 1:(dim-1)];
 
68
  x = permute (x, perm);
 
69
  m = size (x, 1);
 
70
 
 
71
  inc = 1;
 
72
  while (inc*n < m)
 
73
    inc++;
 
74
  endwhile
 
75
  y = fft (x) / m;
 
76
  k = floor (m / 2);
 
77
  sz = size (x);
 
78
  sz(1) = n * inc - m;
 
79
  idx = cell (nd, 1);
 
80
  for i = 2:nd
 
81
    idx{i} = 1:sz(i);
 
82
  endfor
 
83
  idx{1} = 1:k;
 
84
  z = cat (1, y(idx{:}), zeros (sz));
 
85
  idx{1} = k+1:m;
 
86
  z = cat (1, z, y(idx{:}));
 
87
  z = n * ifft (z);
 
88
 
 
89
  if (inc != 1)
 
90
    sz(1) = n;
 
91
    z = inc * reshape (z(1:inc:end), sz);
 
92
  endif
 
93
 
 
94
  z = ipermute (z, perm);
 
95
endfunction
 
96
 
 
97
%!demo
 
98
%! t = 0 : 0.3 : pi; dt = t(2)-t(1);
 
99
%! n = length (t); k = 100;
 
100
%! ti = t(1) + [0 : k-1]*dt*n/k;
 
101
%! y = sin (4*t + 0.3) .* cos (3*t - 0.1);
 
102
%! yp = sin (4*ti + 0.3) .* cos (3*ti - 0.1);
 
103
%! plot (ti, yp, 'g', ti, interp1(t, y, ti, 'spline'), 'b', ...
 
104
%!       ti, interpft (y, k), 'c', t, y, 'r+');
 
105
%! legend ('sin(4t+0.3)cos(3t-0.1','spline','interpft','data');
 
106
 
 
107
%!shared n,y
 
108
%! x = [0:10]'; y = sin(x); n = length (x);
 
109
%!assert (interpft(y, n), y, eps);
 
110
%!assert (interpft(y', n), y', eps);
 
111
%!assert (interpft([y,y],n), [y,y], eps);
 
112
 
 
113
%!error (interpft(y,n,0))
 
114
%!error (interpft(y,[n,n]))