annotate src/DLD-FUNCTIONS/fftn.cc @ 4773:ccfbd6047a54

[project @ 2004-02-16 19:02:32 by jwe]
author jwe
date Mon, 16 Feb 2004 19:02:33 +0000
parents
children 3dfdf6f36854
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
4773
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
1 /*
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
2
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
3 Copyright (C) 2004 David Bateman
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
4
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
6
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
9 Free Software Foundation; either version 2, or (at your option) any
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
10 later version.
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
11
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
15 for more details.
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
16
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
18 along with Octave; see the file COPYING. If not, write to the Free
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
19 Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
20
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
21 */
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
22
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
23 #ifdef HAVE_CONFIG_H
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
24 #include <config.h>
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
25 #endif
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
26
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
27 #include "lo-mappers.h"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
28
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
29 #include "defun-dld.h"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
30 #include "error.h"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
31 #include "gripes.h"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
32 #include "oct-obj.h"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
33 #include "utils.h"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
34
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
35 // This function should be merged with Fifft.
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
36
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
37 #if defined (HAVE_FFTW3)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
38 #define FFTSRC "@sc{Fftw}"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
39 #define WISDOM ", fft_wisdom"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
40 #else
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
41 #define FFTSRC "@sc{Fftpack}"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
42 #define WISDOM
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
43 #endif
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
44
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
45 static octave_value
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
46 do_fftn (const octave_value_list &args, const char *fcn, int type)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
47 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
48 octave_value retval;
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
49
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
50 int nargin = args.length ();
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
51
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
52 if (nargin < 1 || nargin > 2)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
53 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
54 print_usage (fcn);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
55 return retval;
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
56 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
57
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
58 octave_value arg = args(0);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
59 dim_vector dims = arg.dims ();
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
60
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
61 for (int i = 0; i < dims.length (); i++)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
62 if (dims(i) < 0)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
63 return retval;
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
64
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
65 if (nargin > 1)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
66 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
67 Matrix val = args(1).matrix_value ();
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
68 if (val.rows () > val.columns ())
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
69 val.transpose ();
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
70
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
71 if (val.columns () != dims.length () || val.rows () != 1)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
72 error ("%s: second argument must be a vector of length dim", fcn);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
73 else
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
74 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
75 for (int i = 0; i < dims.length (); i++)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
76 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
77 if (xisnan (val(i,0)))
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
78 error ("%s: NaN is invalid as a dimension", fcn);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
79 else if (NINT (val(i,0)) < 0)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
80 error ("%s: all dimension must be greater than zero", fcn);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
81 else
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
82 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
83 dims(i) = NINT(val(i,0));
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
84 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
85 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
86 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
87 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
88
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
89 if (error_state)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
90 return retval;
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
91
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
92 if (dims.all_zero ())
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
93 return octave_value (Matrix ());
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
94
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
95 if (arg.is_real_type ())
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
96 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
97 NDArray nda = arg.array_value ();
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
98
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
99 if (! error_state)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
100 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
101 nda.resize (dims, 0.0);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
102 retval = (type != 0 ? nda.ifourierNd () : nda.fourierNd ());
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
103 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
104 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
105 else if (arg.is_complex_type ())
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
106 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
107 ComplexNDArray cnda = arg.complex_array_value ();
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
108
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
109 if (! error_state)
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
110 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
111 cnda.resize (dims, 0.0);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
112 retval = (type != 0 ? cnda.ifourierNd () : cnda.fourierNd ());
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
113 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
114 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
115 else
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
116 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
117 gripe_wrong_type_arg (fcn, arg);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
118 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
119
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
120 return retval;
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
121 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
122
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
123 DEFUN_DLD (fftn, args, ,
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
124 "-*- texinfo -*-\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
125 @deftypefn {Loadable Function} {} fftn (@var{a}, @var{siz})\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
126 Compute the N dimensional FFT of @var{a} using subroutines from\n"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
127 FFTSRC
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
128 ". The optional vector argument @var{siz} may be used specify the\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
129 dimensions of the array to be used. If an element of @var{siz} is\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
130 smaller than the corresponding dimension, then the dimension is\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
131 truncated prior to performing the FFT. Otherwise if an element\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
132 of @var{siz} is larger than the corresponding dimension @var{a}\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
133 is resized and padded with zeros.\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
134 @end deftypefn\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
135 @seealso {ifftn, fft, fft2"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
136 WISDOM
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
137 "}")
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
138 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
139 return do_fftn(args, "fftn", 0);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
140 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
141
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
142 DEFUN_DLD (ifftn, args, ,
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
143 "-*- texinfo -*-\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
144 @deftypefn {Loadable Function} {} ifftn (@var{a}, @var{siz})\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
145 Compute the invesre N dimensional FFT of @var{a} using subroutines from\n"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
146 FFTSRC
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
147 ". The optional vector argument @var{siz} may be used specify the\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
148 dimensions of the array to be used. If an element of @var{siz} is\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
149 smaller than the corresponding dimension, then the dimension is\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
150 truncated prior to performing the inverse FFT. Otherwise if an element\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
151 of @var{siz} is larger than the corresponding dimension @var{a}\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
152 is resized and padded with zeros.\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
153 @end deftypefn\n\
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
154 @seealso {fftn, ifft, ifft2"
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
155 WISDOM
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
156 "}")
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
157 {
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
158 return do_fftn(args, "ifftn", 1);
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
159 }
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
160
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
161 /*
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
162 ;;; Local Variables: ***
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
163 ;;; mode: C++ ***
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
164 ;;; End: ***
ccfbd6047a54 [project @ 2004-02-16 19:02:32 by jwe]
jwe
parents:
diff changeset
165 */