GNU Octave  4.0.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-2015 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 number of\n\
178 rows and columns of @var{A} to use. If either of these is larger than the\n\
179 size of @var{A}, @var{A} is resized and padded with zeros.\n\
180 \n\
181 If @var{A} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
182 of @var{A} is treated separately.\n\
183 @seealso{ifft2, fft, fftn, fftw}\n\
184 @end deftypefn")
185 {
186  return do_fft2 (args, "fft2", 0);
187 }
188 
189 
190 DEFUN (ifft2, args, ,
191  "-*- texinfo -*-\n\
192 @deftypefn {Built-in Function} {} ifft2 (@var{A})\n\
193 @deftypefnx {Built-in Function} {} ifft2 (@var{A}, @var{m}, @var{n})\n\
194 Compute the inverse two-dimensional discrete Fourier transform of @var{A}\n\
195 using a Fast Fourier Transform (FFT) algorithm.\n\
196 \n\
197 The optional arguments @var{m} and @var{n} may be used specify the number of\n\
198 rows and columns of @var{A} to use. If either of these is larger than the\n\
199 size of @var{A}, @var{A} is resized and padded with zeros.\n\
200 \n\
201 If @var{A} is a multi-dimensional matrix, each two-dimensional sub-matrix\n\
202 of @var{A} is treated separately\n\
203 @seealso{fft2, ifft, ifftn, fftw}\n\
204 @end deftypefn")
205 {
206  return do_fft2 (args, "ifft2", 1);
207 }
208 
209 /*
210 %% Author: David Billinghurst ([email protected])
211 %% Comalco Research and Technology
212 %% 02 May 2000
213 %!test
214 %! M = 16;
215 %! N = 8;
216 %!
217 %! m = 5;
218 %! n = 3;
219 %!
220 %! x = 2*pi*(0:1:M-1)/M;
221 %! y = 2*pi*(0:1:N-1)/N;
222 %! sx = cos (m*x);
223 %! sy = sin (n*y);
224 %! s = kron (sx',sy);
225 %! S = fft2 (s);
226 %! answer = kron (fft (sx)', fft (sy));
227 %! assert (S, answer, 4*M*N*eps);
228 
229 %% Author: David Billinghurst ([email protected])
230 %% Comalco Research and Technology
231 %% 02 May 2000
232 %!test
233 %! M = 12;
234 %! N = 7;
235 %!
236 %! m = 3;
237 %! n = 2;
238 %!
239 %! x = 2*pi*(0:1:M-1)/M;
240 %! y = 2*pi*(0:1:N-1)/N;
241 %!
242 %! sx = cos (m*x);
243 %! sy = cos (n*y);
244 %!
245 %! S = kron (fft (sx)', fft (sy));
246 %! answer = kron (sx', sy);
247 %! s = ifft2 (S);
248 %!
249 %! assert (s, answer, 30*eps);
250 
251 
252 %% Author: David Billinghurst ([email protected])
253 %% Comalco Research and Technology
254 %% 02 May 2000
255 %!test
256 %! M = 16;
257 %! N = 8;
258 %!
259 %! m = 5;
260 %! n = 3;
261 %!
262 %! x = 2*pi*(0:1:M-1)/M;
263 %! y = 2*pi*(0:1:N-1)/N;
264 %! sx = single (cos (m*x));
265 %! sy = single (sin (n*y));
266 %! s = kron (sx', sy);
267 %! S = fft2 (s);
268 %! answer = kron (fft (sx)', fft (sy));
269 %! assert (S, answer, 4*M*N*eps ("single"));
270 
271 %% Author: David Billinghurst ([email protected])
272 %% Comalco Research and Technology
273 %% 02 May 2000
274 %!test
275 %! M = 12;
276 %! N = 7;
277 %!
278 %! m = 3;
279 %! n = 2;
280 %!
281 %! x = single (2*pi*(0:1:M-1)/M);
282 %! y = single (2*pi*(0:1:N-1)/N);
283 %!
284 %! sx = cos (m*x);
285 %! sy = cos (n*y);
286 %!
287 %! S = kron (fft (sx)', fft (sy));
288 %! answer = kron (sx', sy);
289 %! s = ifft2 (S);
290 %!
291 %! assert (s, answer, 30*eps ("single"));
292 */
ComplexNDArray complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:798
bool is_real_type(void) const
Definition: ov.h:651
ComplexNDArray ifourier2d(void) const
Definition: CNDArray.cc:141
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
bool xisnan(double x)
Definition: lo-mappers.cc:144
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
FloatComplexNDArray fourier2d(void) const
Definition: fCNDArray.cc:121
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
static octave_value do_fft2(const octave_value_list &args, const char *fcn, int type)
Definition: fft2.cc:45
ComplexNDArray fourier2d(void) const
Definition: CNDArray.cc:121
FloatComplexNDArray fourier2d(void) const
Definition: fNDArray.cc:119
FloatComplexNDArray ifourier2d(void) const
Definition: fNDArray.cc:139
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:782
FloatComplexNDArray float_complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:802
ComplexNDArray ifourier2d(void) const
Definition: dNDArray.cc:181
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1033
Definition: dMatrix.h:35
dim_vector dims(void) const
Definition: ov.h:470
bool all_zero(void) const
Definition: dim-vector.h:304
double arg(double x)
Definition: lo-mappers.h:37
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:779
ComplexNDArray fourier2d(void) const
Definition: dNDArray.cc:161
FloatComplexNDArray ifourier2d(void) const
Definition: fCNDArray.cc:141
bool is_single_type(void) const
Definition: ov.h:611
octave_idx_type NINTbig(double x)
Definition: lo-mappers.cc:635
int length(void) const
Definition: dim-vector.h:281
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))