GNU Octave  4.2.1
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-2017 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 #if defined (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 "errwarn.h"
33 #include "ovl.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  int nargin = args.length ();
46 
47  if (nargin < 1 || nargin > 3)
48  print_usage ();
49 
51  octave_value arg = args(0);
52  octave_idx_type n_points = -1;
53  dim_vector dims = arg.dims ();
54  int ndims = dims.ndims ();
55  int dim = -1;
56 
57  if (nargin > 1)
58  {
59  if (! args(1).is_empty ())
60  {
61  double dval = args(1).double_value ();
62  if (octave::math::isnan (dval))
63  error ("%s: number of points (N) cannot be NaN", fcn);
64 
65  n_points = octave::math::nint_big (dval);
66  if (n_points < 0)
67  error ("%s: number of points (N) must be greater than zero", fcn);
68  }
69  }
70 
71  if (nargin > 2)
72  {
73  double dval = args(2).double_value ();
74  if (octave::math::isnan (dval))
75  error ("%s: DIM cannot be NaN", fcn);
76  else if (dval < 1 || dval > ndims)
77  error ("%s: DIM must be a valid dimension along which to perform FFT",
78  fcn);
79  else
80  // to be safe, cast it back to int since dim is an int
81  dim = octave::math::nint (dval) - 1;
82  }
83 
84  // FIXME: This seems strange and unnecessary (10/21/16).
85  // How would you ever arrive at an octave_value object without correct dims?
86  // We certainly don't make this check every other place in Octave.
87  for (octave_idx_type i = 0; i < ndims; i++)
88  if (dims(i) < 0)
89  return retval;
90 
91  if (dim < 0)
92  {
93  dim = dims.first_non_singleton ();
94 
95  // And if the first argument is scalar?
96  if (dim == ndims)
97  dim = 1;
98  }
99 
100  if (n_points < 0)
101  n_points = dims(dim);
102  else
103  dims(dim) = n_points;
104 
105  if (n_points == 0 || dims.any_zero ())
106  {
107  if (arg.is_single_type ())
108  return octave_value (FloatNDArray (dims));
109  else
110  return octave_value (NDArray (dims));
111  }
112 
113  if (n_points == 1)
114  {
115  octave_value_list idx (ndims);
116  for (octave_idx_type i = 0; i < ndims; i++)
117  idx(i) = idx_vector::colon;
118  idx(dim) = idx_vector (static_cast<octave_idx_type> (0));
119 
120  return arg.do_index_op (idx);
121  }
122 
123  if (arg.is_single_type ())
124  {
125  if (arg.is_real_type ())
126  {
127  FloatNDArray nda = arg.float_array_value ();
128 
129  nda.resize (dims, 0.0);
130  retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
131  }
132  else
133  {
135 
136  cnda.resize (dims, 0.0);
137  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
138  }
139  }
140  else
141  {
142  if (arg.is_real_type ())
143  {
144  NDArray nda = arg.array_value ();
145 
146  nda.resize (dims, 0.0);
147  retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim));
148  }
149  else if (arg.is_complex_type ())
150  {
151  ComplexNDArray cnda = arg.complex_array_value ();
152 
153  cnda.resize (dims, 0.0);
154  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
155  }
156  else
157  err_wrong_type_arg (fcn, arg);
158  }
159 
160  return retval;
161 }
162 
163 /*
164 %!assert (fft ([]), [])
165 %!assert (fft (zeros (10,0)), zeros (10,0))
166 %!assert (fft (zeros (0,10)), zeros (0,10))
167 %!assert (fft (0), 0)
168 %!assert (fft (1), 1)
169 %!assert (fft (ones (2,2)), [2,2; 0,0])
170 %!assert (fft (eye (2,2)), [1,1; 1,-1])
171 
172 %!assert (fft (single ([])), single ([]))
173 %!assert (fft (zeros (10,0,"single")), zeros (10,0,"single"))
174 %!assert (fft (zeros (0,10,"single")), zeros (0,10,"single"))
175 %!assert (fft (single (0)), single (0))
176 %!assert (fft (single (1)), single (1))
177 %!assert (fft (ones (2,2,"single")), single ([2,2; 0,0]))
178 %!assert (fft (eye (2,2,"single")), single ([1,1; 1,-1]))
179 
180 %!error (fft ())
181 */
182 
183 
184 DEFUN (fft, args, ,
185  doc: /* -*- texinfo -*-
186 @deftypefn {} {} fft (@var{x})
187 @deftypefnx {} {} fft (@var{x}, @var{n})
188 @deftypefnx {} {} fft (@var{x}, @var{n}, @var{dim})
189 Compute the discrete Fourier transform of @var{A} using
190 a Fast Fourier Transform (FFT) algorithm.
191 
192 The FFT is calculated along the first non-singleton dimension of the
193 array. Thus if @var{x} is a matrix, @code{fft (@var{x})} computes the
194 FFT for each column of @var{x}.
195 
196 If called with two arguments, @var{n} is expected to be an integer
197 specifying the number of elements of @var{x} to use, or an empty
198 matrix to specify that its value should be ignored. If @var{n} is
199 larger than the dimension along which the FFT is calculated, then
200 @var{x} is resized and padded with zeros. Otherwise, if @var{n} is
201 smaller than the dimension along which the FFT is calculated, then
202 @var{x} is truncated.
203 
204 If called with three arguments, @var{dim} is an integer specifying the
205 dimension of the matrix along which the FFT is performed
206 @seealso{ifft, fft2, fftn, fftw}
207 @end deftypefn */)
208 {
209  return do_fft (args, "fft", 0);
210 }
211 
212 
213 DEFUN (ifft, args, ,
214  doc: /* -*- texinfo -*-
215 @deftypefn {} {} ifft (@var{x})
216 @deftypefnx {} {} ifft (@var{x}, @var{n})
217 @deftypefnx {} {} ifft (@var{x}, @var{n}, @var{dim})
218 Compute the inverse discrete Fourier transform of @var{A}
219 using a Fast Fourier Transform (FFT) algorithm.
220 
221 The inverse FFT is calculated along the first non-singleton dimension
222 of the array. Thus if @var{x} is a matrix, @code{fft (@var{x})} computes
223 the inverse FFT for each column of @var{x}.
224 
225 If called with two arguments, @var{n} is expected to be an integer
226 specifying the number of elements of @var{x} to use, or an empty
227 matrix to specify that its value should be ignored. If @var{n} is
228 larger than the dimension along which the inverse FFT is calculated, then
229 @var{x} is resized and padded with zeros. Otherwise, if @var{n} is
230 smaller than the dimension along which the inverse FFT is calculated,
231 then @var{x} is truncated.
232 
233 If called with three arguments, @var{dim} is an integer specifying the
234 dimension of the matrix along which the inverse FFT is performed
235 @seealso{fft, ifft2, ifftn, fftw}
236 @end deftypefn */)
237 {
238  return do_fft (args, "ifft", 1);
239 }
240 
241 /*
242 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
243 ## Comalco Research and Technology
244 ## 02 May 2000
245 %!test
246 %! N = 64;
247 %! n = 4;
248 %! t = 2*pi*(0:1:N-1)/N;
249 %! s = cos (n*t);
250 %! S = fft (s);
251 %!
252 %! answer = zeros (size (t));
253 %! answer(n+1) = N/2;
254 %! answer(N-n+1) = N/2;
255 %!
256 %! assert (S, answer, 4*N*eps);
257 
258 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
259 ## Comalco Research and Technology
260 ## 02 May 2000
261 %!test
262 %! N = 64;
263 %! n = 7;
264 %! t = 2*pi*(0:1:N-1)/N;
265 %! s = cos (n*t);
266 %!
267 %! S = zeros (size (t));
268 %! S(n+1) = N/2;
269 %! S(N-n+1) = N/2;
270 %!
271 %! assert (ifft (S), s, 4*N*eps);
272 
273 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
274 ## Comalco Research and Technology
275 ## 02 May 2000
276 %!test
277 %! N = 64;
278 %! n = 4;
279 %! t = single (2*pi*(0:1:N-1)/N);
280 %! s = cos (n*t);
281 %! S = fft (s);
282 %!
283 %! answer = zeros (size (t), "single");
284 %! answer(n+1) = N/2;
285 %! answer(N-n+1) = N/2;
286 %!
287 %! assert (S, answer, 4*N*eps ("single"));
288 
289 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
290 ## Comalco Research and Technology
291 ## 02 May 2000
292 %!test
293 %! N = 64;
294 %! n = 7;
295 %! t = 2*pi*(0:1:N-1)/N;
296 %! s = cos (n*t);
297 %!
298 %! S = zeros (size (t), "single");
299 %! S(n+1) = N/2;
300 %! S(N-n+1) = N/2;
301 %!
302 %! assert (ifft (S), s, 4*N*eps ("single"));
303 */
FloatComplexNDArray ifourier(int dim=1) const
Definition: fCNDArray.cc:89
FloatComplexNDArray fourier(int dim=1) const
Definition: fNDArray.cc:58
ComplexNDArray complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:812
bool is_real_type(void) const
Definition: ov.h:667
static const idx_vector colon
Definition: idx-vector.h:482
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
octave_idx_type length(void) const
Definition: ovl.h:96
bool isnan(double x)
Definition: lo-mappers.cc:347
ComplexNDArray fourier(int dim=1) const
Definition: CNDArray.cc:58
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
FloatComplexNDArray ifourier(int dim=1) const
Definition: fNDArray.cc:89
int first_non_singleton(int def=0) const
Definition: dim-vector.h:463
octave_value arg
Definition: pr-output.cc:3440
FloatComplexNDArray fourier(int dim=1) const
Definition: fCNDArray.cc:58
octave_function * fcn
Definition: ov-class.cc:1743
static octave_value do_fft(const octave_value_list &args, const char *fcn, int type)
Definition: fft.cc:43
ComplexNDArray ifourier(int dim=1) const
Definition: CNDArray.cc:89
JNIEnv void * args
Definition: ov-java.cc:67
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:796
ComplexNDArray ifourier(int dim=1) const
Definition: dNDArray.cc:131
int nargin
Definition: graphics.cc:10115
FloatComplexNDArray float_complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:816
bool is_complex_type(void) const
Definition: ov.h:670
void resize(const dim_vector &dv, const T &rfv)
Definition: Array.cc:1028
octave_value retval
Definition: data.cc:6294
idx type
Definition: ov.cc:3129
bool any_zero(void) const
Definition: dim-vector.h:359
the exceeded dimensions are set to if fewer subscripts than dimensions are the exceeding dimensions are merged into the final requested dimension For consider the following dims
Definition: sub2ind.cc:255
dim_vector dims(void) const
Definition: ov.h:486
void err_wrong_type_arg(const char *name, const char *s)
Definition: errwarn.cc:156
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:793
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave_idx_type nint_big(double x)
Definition: lo-mappers.cc:409
ComplexNDArray fourier(int dim=1) const
Definition: dNDArray.cc:100
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:301
bool is_single_type(void) const
Definition: ov.h:627
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
int nint(double x)
Definition: lo-mappers.cc:433
octave_value do_index_op(const octave_value_list &idx, bool resize_ok=false)
Definition: ov.h:454