Mercurial > hg > octave-nkf
annotate scripts/general/interpft.m @ 9665:1dba57e9d08d
use blas_trans_type for xgemm
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Sat, 26 Sep 2009 10:41:07 +0200 |
| parents | 1bf0ce0930be |
| children | be55736a0783 |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2001, 2006, 2007, 2008, 2009 Paul Kienzle |
| 5837 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 5837 | 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 | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 5837 | 18 |
| 19 ## -*- texinfo -*- | |
| 20 ## @deftypefn {Function File} {} interpft (@var{x}, @var{n}) | |
| 21 ## @deftypefnx {Function File} {} interpft (@var{x}, @var{n}, @var{dim}) | |
| 22 ## | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
23 ## Fourier interpolation. If @var{x} is a vector, then @var{x} is |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
24 ## resampled with @var{n} points. The data in @var{x} is assumed to be |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
25 ## equispaced. If @var{x} is an array, then operate along each column of |
|
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
26 ## the array separately. If @var{dim} is specified, then interpolate |
| 5837 | 27 ## along the dimension @var{dim}. |
| 28 ## | |
| 29 ## @code{interpft} assumes that the interpolated function is periodic, | |
| 7001 | 30 ## and so assumptions are made about the end points of the interpolation. |
| 5837 | 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) | |
| 5838 | 50 if (isvector (x) && size (x, 1) == 1) |
| 5837 | 51 dim = 2; |
| 52 else | |
| 53 dim = 1; | |
| 54 endif | |
| 55 endif | |
| 56 | |
| 5838 | 57 if (! isscalar (n)) |
| 5837 | 58 error ("interpft: n must be an integer scalar"); |
| 59 endif | |
| 60 | |
| 5838 | 61 nd = ndims (x); |
| 5837 | 62 |
| 63 if (dim < 1 || dim > nd) | |
| 64 error ("interpft: integrating over invalid dimension"); | |
| 65 endif | |
| 66 | |
| 5838 | 67 perm = [dim:nd, 1:(dim-1)]; |
| 68 x = permute (x, perm); | |
| 5837 | 69 m = size (x, 1); |
| 70 | |
| 71 inc = 1; | |
| 5838 | 72 while (inc*n < m) |
| 5837 | 73 inc++; |
| 74 endwhile | |
| 75 y = fft (x) / m; | |
| 76 k = floor (m / 2); | |
| 5838 | 77 sz = size (x); |
| 5837 | 78 sz(1) = n * inc - m; |
| 5838 | 79 idx = cell (nd, 1); |
| 5837 | 80 for i = 2:nd |
| 81 idx{i} = 1:sz(i); | |
| 82 endfor | |
| 83 idx{1} = 1:k; | |
| 5838 | 84 z = cat (1, y(idx{:}), zeros (sz)); |
| 5837 | 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; | |
| 5838 | 91 z = inc * reshape (z(1:inc:end), sz); |
| 5837 | 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); | |
| 6702 | 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'); | |
| 5837 | 106 |
| 107 %!shared n,y | |
| 108 %! x = [0:10]'; y = sin(x); n = length (x); | |
|
8656
e3041433a57e
increase tolerance in interpft tests
Rob Mahurin <rob@chimera.(none)>
parents:
8138
diff
changeset
|
109 %!assert (interpft(y, n), y, 20*eps); |
|
e3041433a57e
increase tolerance in interpft tests
Rob Mahurin <rob@chimera.(none)>
parents:
8138
diff
changeset
|
110 %!assert (interpft(y', n), y', 20*eps); |
|
e3041433a57e
increase tolerance in interpft tests
Rob Mahurin <rob@chimera.(none)>
parents:
8138
diff
changeset
|
111 %!assert (interpft([y,y],n), [y,y], 20*eps); |
| 5837 | 112 |
| 113 %!error (interpft(y,n,0)) | |
| 114 %!error (interpft(y,[n,n])) |