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
fft.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 #if defined (HAVE_FFTW)
37 #define FFTSRC "@sc{fftw}"
38 #else
39 #define FFTSRC "@sc{fftpack}"
40 #endif
41 
42 static octave_value
43 do_fft (const octave_value_list &args, const char *fcn, int type)
44 {
45  octave_value retval;
46 
47  int nargin = args.length ();
48 
49  if (nargin < 1 || nargin > 3)
50  {
51  print_usage ();
52  return retval;
53  }
54 
55  octave_value arg = args(0);
56  dim_vector dims = arg.dims ();
57  octave_idx_type n_points = -1;
58  int dim = -1;
59 
60  if (nargin > 1)
61  {
62  if (! args(1).is_empty ())
63  {
64  double dval = args(1).double_value ();
65  if (xisnan (dval))
66  error ("%s: number of points (N) cannot be NaN", fcn);
67  else
68  {
69  n_points = NINTbig (dval);
70  if (n_points < 0)
71  error ("%s: number of points (N) must be greater than zero",
72  fcn);
73  }
74  }
75  }
76 
77  if (error_state)
78  return retval;
79 
80  if (nargin > 2)
81  {
82  double dval = args(2).double_value ();
83  if (xisnan (dval))
84  error ("%s: DIM cannot be NaN", fcn);
85  else if (dval < 1 || dval > dims.length ())
86  error ("%s: DIM must be a valid dimension along which to perform FFT",
87  fcn);
88  else
89  // to be safe, cast it back to int since dim is an int
90  dim = NINT (dval) - 1;
91  }
92 
93  if (error_state)
94  return retval;
95 
96  for (octave_idx_type i = 0; i < dims.length (); i++)
97  if (dims(i) < 0)
98  return retval;
99 
100  if (dim < 0)
101  {
102  for (octave_idx_type i = 0; i < dims.length (); i++)
103  if (dims(i) > 1)
104  {
105  dim = i;
106  break;
107  }
108 
109  // And if the first argument is scalar?
110  if (dim < 0)
111  dim = 1;
112  }
113 
114  if (n_points < 0)
115  n_points = dims (dim);
116  else
117  dims (dim) = n_points;
118 
119  if (dims.any_zero () || n_points == 0)
120  {
121  if (arg.is_single_type ())
122  return octave_value (FloatNDArray (dims));
123  else
124  return octave_value (NDArray (dims));
125  }
126 
127  if (arg.is_single_type ())
128  {
129  if (arg.is_real_type ())
130  {
131  FloatNDArray nda = arg.float_array_value ();
132 
133  if (! error_state)
134  {
135  nda.resize (dims, 0.0);
136  retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
137  }
138  }
139  else
140  {
142 
143  if (! error_state)
144  {
145  cnda.resize (dims, 0.0);
146  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
147  }
148  }
149  }
150  else
151  {
152  if (arg.is_real_type ())
153  {
154  NDArray nda = arg.array_value ();
155 
156  if (! error_state)
157  {
158  nda.resize (dims, 0.0);
159  retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
160  }
161  }
162  else if (arg.is_complex_type ())
163  {
164  ComplexNDArray cnda = arg.complex_array_value ();
165 
166  if (! error_state)
167  {
168  cnda.resize (dims, 0.0);
169  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
170  }
171  }
172  else
173  {
174  gripe_wrong_type_arg (fcn, arg);
175  }
176  }
177 
178  return retval;
179 }
180 
181 /*
182 %!assert (fft ([]), [])
183 %!assert (fft (zeros (10,0)), zeros (10,0))
184 %!assert (fft (zeros (0,10)), zeros (0,10))
185 %!assert (fft (0), 0)
186 %!assert (fft (1), 1)
187 %!assert (fft (ones (2,2)), [2,2; 0,0])
188 %!assert (fft (eye (2,2)), [1,1; 1,-1])
189 
190 %!assert (fft (single ([])), single ([]))
191 %!assert (fft (zeros (10,0,"single")), zeros (10,0,"single"))
192 %!assert (fft (zeros (0,10,"single")), zeros (0,10,"single"))
193 %!assert (fft (single (0)), single (0))
194 %!assert (fft (single (1)), single (1))
195 %!assert (fft (ones (2,2,"single")), single ([2,2; 0,0]))
196 %!assert (fft (eye (2,2,"single")), single ([1,1; 1,-1]))
197 
198 %!error (fft ())
199 */
200 
201 
202 DEFUN (fft, args, ,
203  "-*- texinfo -*-\n\
204 @deftypefn {Built-in Function} {} fft (@var{x})\n\
205 @deftypefnx {Built-in Function} {} fft (@var{x}, @var{n})\n\
206 @deftypefnx {Built-in Function} {} fft (@var{x}, @var{n}, @var{dim})\n\
207 Compute the discrete Fourier transform of @var{A} using\n\
208 a Fast Fourier Transform (FFT) algorithm.\n\
209 \n\
210 The FFT is calculated along the first non-singleton dimension of the\n\
211 array. Thus if @var{x} is a matrix, @code{fft (@var{x})} computes the\n\
212 FFT for each column of @var{x}.\n\
213 \n\
214 If called with two arguments, @var{n} is expected to be an integer\n\
215 specifying the number of elements of @var{x} to use, or an empty\n\
216 matrix to specify that its value should be ignored. If @var{n} is\n\
217 larger than the dimension along which the FFT is calculated, then\n\
218 @var{x} is resized and padded with zeros. Otherwise, if @var{n} is\n\
219 smaller than the dimension along which the FFT is calculated, then\n\
220 @var{x} is truncated.\n\
221 \n\
222 If called with three arguments, @var{dim} is an integer specifying the\n\
223 dimension of the matrix along which the FFT is performed\n\
224 @seealso{ifft, fft2, fftn, fftw}\n\
225 @end deftypefn")
226 {
227  return do_fft (args, "fft", 0);
228 }
229 
230 
231 DEFUN (ifft, args, ,
232  "-*- texinfo -*-\n\
233 @deftypefn {Built-in Function} {} ifft (@var{x})\n\
234 @deftypefnx {Built-in Function} {} ifft (@var{x}, @var{n})\n\
235 @deftypefnx {Built-in Function} {} ifft (@var{x}, @var{n}, @var{dim})\n\
236 Compute the inverse discrete Fourier transform of @var{A}\n\
237 using a Fast Fourier Transform (FFT) algorithm.\n\
238 \n\
239 The inverse FFT is calculated along the first non-singleton dimension\n\
240 of the array. Thus if @var{x} is a matrix, @code{fft (@var{x})} computes\n\
241 the inverse FFT for each column of @var{x}.\n\
242 \n\
243 If called with two arguments, @var{n} is expected to be an integer\n\
244 specifying the number of elements of @var{x} to use, or an empty\n\
245 matrix to specify that its value should be ignored. If @var{n} is\n\
246 larger than the dimension along which the inverse FFT is calculated, then\n\
247 @var{x} is resized and padded with zeros. Otherwise, if @var{n} is\n\
248 smaller than the dimension along which the inverse FFT is calculated,\n\
249 then @var{x} is truncated.\n\
250 \n\
251 If called with three arguments, @var{dim} is an integer specifying the\n\
252 dimension of the matrix along which the inverse FFT is performed\n\
253 @seealso{fft, ifft2, ifftn, fftw}\n\
254 @end deftypefn")
255 {
256  return do_fft (args, "ifft", 1);
257 }
258 
259 /*
260 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
261 %% Comalco Research and Technology
262 %% 02 May 2000
263 %!test
264 %! N = 64;
265 %! n = 4;
266 %! t = 2*pi*(0:1:N-1)/N;
267 %! s = cos (n*t);
268 %! S = fft (s);
269 %!
270 %! answer = zeros (size (t));
271 %! answer(n+1) = N/2;
272 %! answer(N-n+1) = N/2;
273 %!
274 %! assert (S, answer, 4*N*eps);
275 
276 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
277 %% Comalco Research and Technology
278 %% 02 May 2000
279 %!test
280 %! N = 64;
281 %! n = 7;
282 %! t = 2*pi*(0:1:N-1)/N;
283 %! s = cos (n*t);
284 %!
285 %! S = zeros (size (t));
286 %! S(n+1) = N/2;
287 %! S(N-n+1) = N/2;
288 %!
289 %! assert (ifft (S), s, 4*N*eps);
290 
291 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
292 %% Comalco Research and Technology
293 %% 02 May 2000
294 %!test
295 %! N = 64;
296 %! n = 4;
297 %! t = single (2*pi*(0:1:N-1)/N);
298 %! s = cos (n*t);
299 %! S = fft (s);
300 %!
301 %! answer = zeros (size (t), "single");
302 %! answer(n+1) = N/2;
303 %! answer(N-n+1) = N/2;
304 %!
305 %! assert (S, answer, 4*N*eps ("single"));
306 
307 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
308 %% Comalco Research and Technology
309 %% 02 May 2000
310 %!test
311 %! N = 64;
312 %! n = 7;
313 %! t = 2*pi*(0:1:N-1)/N;
314 %! s = cos (n*t);
315 %!
316 %! S = zeros (size (t), "single");
317 %! S(n+1) = N/2;
318 %! S(N-n+1) = N/2;
319 %!
320 %! assert (ifft (S), s, 4*N*eps ("single"));
321 */