Mercurial > hg > octave-nkf

comparison src/corefcn/fft2.cc @ 15039:e753177cde93

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory * __contourc__.cc, __dispatch__.cc, __lin_interpn__.cc, __pchip_deriv__.cc, __qp__.cc, balance.cc, besselj.cc, betainc.cc, bsxfun.cc, cellfun.cc, colloc.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, det.cc, dlmread.cc, dot.cc, eig.cc, fft.cc, fft2.cc, fftn.cc, filter.cc, find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc, hess.cc, hex2num.cc, inv.cc, kron.cc, lookup.cc, lsode.cc, lu.cc, luinc.cc, matrix_type.cc, max.cc, md5sum.cc, mgorth.cc, nproc.cc, pinv.cc, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, spparms.cc, sqrtm.cc, str2double.cc, strfind.cc, sub2ind.cc, svd.cc, syl.cc, time.cc, tril.cc, typecast.cc: Move functions from DLD-FUNCTIONS/ to corefcn/ directory. Include "defun.h", not "defun-dld.h". Change docstring to refer to these as "Built-in Functions". * build-aux/mk-opts.pl: Generate options code with '#include "defun.h"'. Change option docstrings to refer to these as "Built-in Functions". * corefcn/module.mk: List of functions to build in corefcn/ dir. * DLD-FUNCTIONS/config-module.awk: Update to new build system. * DLD-FUNCTIONS/module-files: Remove functions which are now in corefcn/ directory. * src/Makefile.am: Update to build "convenience library" in corefcn/. Octave program now links against all other libraries + corefcn libary. * src/find-defun-files.sh: Strip $srcdir from filename. * src/link-deps.mk: Add REGEX and FFTW link dependencies for liboctinterp. * type.m, which.m: Change failing tests to use 'amd', still a dynamic function, rather than 'dot', which isn't.
author Rik <rik@octave.org>
date Fri, 27 Jul 2012 15:35:00 -0700
parents src/DLD-FUNCTIONS/fft2.cc@60e5cf354d80
children
comparison
equal deleted inserted replaced
15038:ab18578c2ade 15039:e753177cde93
1 /*
2
3 Copyright (C) 1997-2012 David Bateman
4 Copyright (C) 1996-1997 John W. Eaton
5
6 This file is part of Octave.
7
8 Octave is free software; you can redistribute it and/or modify it
9 under the terms of the GNU General Public License as published by the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21
22 */
23
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27
28 #include "lo-mappers.h"
29
30 #include "defun.h"
31 #include "error.h"
32 #include "gripes.h"
33 #include "oct-obj.h"
34 #include "utils.h"
35
36 // This function should be merged with Fifft.
37
38 #if defined (HAVE_FFTW)
39 #define FFTSRC "@sc{fftw}"
40 #else
41 #define FFTSRC "@sc{fftpack}"
42 #endif
43
44 static octave_value
45 do_fft2 (const octave_value_list &args, const char *fcn, int type)
46 {
47 octave_value retval;
48
49 int nargin = args.length ();
50
51 if (nargin < 1 || nargin > 3)
52 {
53 print_usage ();
54 return retval;
55 }
56
57 octave_value arg = args(0);
58 dim_vector dims = arg.dims ();
59 octave_idx_type n_rows = -1;
60
61 if (nargin > 1)
62 {
63 double dval = args(1).double_value ();
64 if (xisnan (dval))
65 error ("%s: number of rows (N) cannot be NaN", fcn);
66 else
67 {
68 n_rows = NINTbig (dval);
69 if (n_rows < 0)
70 error ("%s: number of rows (N) must be greater than zero", fcn);
71 }
72 }
73
74 if (error_state)
75 return retval;
76
77 octave_idx_type n_cols = -1;
78 if (nargin > 2)
79 {
80 double dval = args(2).double_value ();
81 if (xisnan (dval))
82 error ("%s: number of columns (M) cannot be NaN", fcn);
83 else
84 {
85 n_cols = NINTbig (dval);
86 if (n_cols < 0)
87 error ("%s: number of columns (M) must be greater than zero", fcn);
88 }
89 }
90
91 if (error_state)
92 return retval;
93
94 for (int i = 0; i < dims.length (); i++)
95 if (dims(i) < 0)
96 return retval;
97
98 if (n_rows < 0)
99 n_rows = dims (0);
100 else
101 dims (0) = n_rows;
102
103 if (n_cols < 0)
104 n_cols = dims (1);
105 else
106 dims (1) = n_cols;
107
108 if (dims.all_zero () || n_rows == 0 || n_cols == 0)
109 {
110 if (arg.is_single_type ())
111 return octave_value (FloatMatrix ());
112 else
113 return octave_value (Matrix ());
114 }
115
116 if (arg.is_single_type ())
117 {
118 if (arg.is_real_type ())
119 {
120 FloatNDArray nda = arg.float_array_value ();
121
122 if (! error_state)
123 {
124 nda.resize (dims, 0.0);
125 retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
126 }
127 }
128 else
129 {
130 FloatComplexNDArray cnda = arg.float_complex_array_value ();
131
132 if (! error_state)
133 {
134 cnda.resize (dims, 0.0);
135 retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
136 }
137 }
138 }
139 else
140 {
141 if (arg.is_real_type ())
142 {
143 NDArray nda = arg.array_value ();
144
145 if (! error_state)
146 {
147 nda.resize (dims, 0.0);
148 retval = (type != 0 ? nda.ifourier2d () : nda.fourier2d ());
149 }
150 }
151 else if (arg.is_complex_type ())
152 {
153 ComplexNDArray cnda = arg.complex_array_value ();
154
155 if (! error_state)
156 {
157 cnda.resize (dims, 0.0);
158 retval = (type != 0 ? cnda.ifourier2d () : cnda.fourier2d ());
159 }
160 }
161 else
162 {
163 gripe_wrong_type_arg (fcn, arg);
164 }
165 }
166
167 return retval;
168 }
169
170 DEFUN (fft2, args, ,
171 "-*- texinfo -*-\n\
172 @deftypefn {Built-in Function} {} fft2 (@var{A})\n\
173 @deftypefnx {Built-in Function} {} fft2 (@var{A}, @var{m}, @var{n})\n\
174 Compute the two-dimensional discrete Fourier transform of @var{A} using\n\
175 a Fast Fourier Transform (FFT) algorithm.\n\
176 \n\
177 The optional arguments @var{m} and @var{n} may be used specify the\n\
178 number of rows and columns of @var{A} to use. If either of these is\n\
179 larger than the size of @var{A}, @var{A} is resized and padded with\n\
180 zeros.\n\
181 \n\
182 If @var{A} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
183 of @var{A} is treated separately.\n\
184 @seealso {ifft2, fft, fftn, fftw}\n\
185 @end deftypefn")
186 {
187 return do_fft2 (args, "fft2", 0);
188 }
189
190
191 DEFUN (ifft2, args, ,
192 "-*- texinfo -*-\n\
193 @deftypefn {Built-in Function} {} ifft2 (@var{A})\n\
194 @deftypefnx {Built-in Function} {} ifft2 (@var{A}, @var{m}, @var{n})\n\
195 Compute the inverse two-dimensional discrete Fourier transform of @var{A}\n\
196 using a Fast Fourier Transform (FFT) algorithm.\n\
197 \n\
198 The optional arguments @var{m} and @var{n} may be used specify the\n\
199 number of rows and columns of @var{A} to use. If either of these is\n\
200 larger than the size of @var{A}, @var{A} is resized and padded with\n\
201 zeros.\n\
202 \n\
203 If @var{A} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
204 of @var{A} is treated separately\n\
205 @seealso {fft2, ifft, ifftn, fftw}\n\
206 @end deftypefn")
207 {
208 return do_fft2 (args, "ifft2", 1);
209 }
210
211 /*
212 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
213 %% Comalco Research and Technology
214 %% 02 May 2000
215 %!test
216 %! M = 16;
217 %! N = 8;
218 %!
219 %! m = 5;
220 %! n = 3;
221 %!
222 %! x = 2*pi*(0:1:M-1)/M;
223 %! y = 2*pi*(0:1:N-1)/N;
224 %! sx = cos (m*x);
225 %! sy = sin (n*y);
226 %! s = kron (sx',sy);
227 %! S = fft2 (s);
228 %! answer = kron (fft (sx)', fft (sy));
229 %! assert (S, answer, 4*M*N*eps);
230
231 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
232 %% Comalco Research and Technology
233 %% 02 May 2000
234 %!test
235 %! M = 12;
236 %! N = 7;
237 %!
238 %! m = 3;
239 %! n = 2;
240 %!
241 %! x = 2*pi*(0:1:M-1)/M;
242 %! y = 2*pi*(0:1:N-1)/N;
243 %!
244 %! sx = cos (m*x);
245 %! sy = cos (n*y);
246 %!
247 %! S = kron (fft (sx)', fft (sy));
248 %! answer = kron (sx', sy);
249 %! s = ifft2 (S);
250 %!
251 %! assert (s, answer, 30*eps);
252
253
254 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
255 %% Comalco Research and Technology
256 %% 02 May 2000
257 %!test
258 %! M = 16;
259 %! N = 8;
260 %!
261 %! m = 5;
262 %! n = 3;
263 %!
264 %! x = 2*pi*(0:1:M-1)/M;
265 %! y = 2*pi*(0:1:N-1)/N;
266 %! sx = single (cos (m*x));
267 %! sy = single (sin (n*y));
268 %! s = kron (sx', sy);
269 %! S = fft2 (s);
270 %! answer = kron (fft (sx)', fft (sy));
271 %! assert (S, answer, 4*M*N*eps ("single"));
272
273 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
274 %% Comalco Research and Technology
275 %% 02 May 2000
276 %!test
277 %! M = 12;
278 %! N = 7;
279 %!
280 %! m = 3;
281 %! n = 2;
282 %!
283 %! x = single (2*pi*(0:1:M-1)/M);
284 %! y = single (2*pi*(0:1:N-1)/N);
285 %!
286 %! sx = cos (m*x);
287 %! sy = cos (n*y);
288 %!
289 %! S = kron (fft (sx)', fft (sy));
290 %! answer = kron (sx', sy);
291 %! s = ifft2 (S);
292 %!
293 %! assert (s, answer, 30*eps ("single"));
294 */