GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
fft2.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1997-2013 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  {
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 ([email protected])
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 ([email protected])
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 ([email protected])
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 ([email protected])
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 */