GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
fft.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1997-2018 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
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License 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 <https://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 
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).isempty ())
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.isreal ())
126  {
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.isreal ())
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.iscomplex ())
150  {
152 
153  cnda.resize (dims, 0.0);
154  retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim));
155  }
156  else
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{x} 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{x}
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 */
ComplexNDArray fourier(int dim=1) const
Definition: dNDArray.cc:98
ComplexNDArray ifourier(int dim=1) const
Definition: dNDArray.cc:129
static const idx_vector colon
Definition: idx-vector.h:498
OCTINTERP_API void print_usage(void)
Definition: defun.cc:54
bool isnan(bool)
Definition: lo-mappers.h:187
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:53
void error(const char *fmt,...)
Definition: error.cc:578
FloatComplexNDArray fourier(int dim=1) const
Definition: fNDArray.cc:56
FloatNDArray float_array_value(bool frc_str_conv=false) const
Definition: ov.h:843
FloatComplexNDArray ifourier(int dim=1) const
Definition: fCNDArray.cc:87
octave_value arg
Definition: pr-output.cc:3244
octave_function * fcn
Definition: ov-class.cc:1754
ComplexNDArray fourier(int dim=1) const
Definition: CNDArray.cc:56
static octave_value do_fft(const octave_value_list &args, const char *fcn, int type)
Definition: fft.cc:43
bool is_single_type(void) const
Definition: ov.h:651
dim_vector dims(void) const
Definition: ov.h:469
FloatComplexNDArray fourier(int dim=1) const
Definition: fCNDArray.cc:56
void resize(const dim_vector &dv, const T &rfv)
Resizing (with fill).
Definition: Array.cc:1010
octave_value retval
Definition: data.cc:6246
idx type
Definition: ov.cc:3114
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
FloatComplexNDArray float_complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:863
void err_wrong_type_arg(const char *name, const char *s)
Definition: errwarn.cc:162
return octave_value(v1.char_array_value() . concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string()) ? '\'' :'"'))
bool isreal(void) const
Definition: ov.h:703
octave_idx_type nint_big(double x)
Definition: lo-mappers.cc:182
octave_idx_type length(void) const
Definition: ovl.h:96
args.length() nargin
Definition: file-io.cc:589
bool iscomplex(void) const
Definition: ov.h:710
for i
Definition: data.cc:5264
ComplexNDArray ifourier(int dim=1) const
Definition: CNDArray.cc:87
ComplexNDArray complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:859
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
FloatComplexNDArray ifourier(int dim=1) const
Definition: fNDArray.cc:87
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:840
int nint(double x)
Definition: lo-mappers.cc:206
octave_value do_index_op(const octave_value_list &idx, bool resize_ok=false)
Definition: ov.h:444