Mercurial > octave-nkf
view src/DLD-FUNCTIONS/fft.cc @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | c827f5673321 |
children | 87865ed7405f |
line wrap: on
line source
/* Copyright (C) 1997, 1999, 2002, 2004, 2005, 2006, 2007 David Bateman Copyright (C) 1996, 1997 John W. Eaton This file is part of Octave. Octave is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. */ #ifdef HAVE_CONFIG_H #include <config.h> #endif #include "lo-mappers.h" #include "defun-dld.h" #include "error.h" #include "gripes.h" #include "oct-obj.h" #include "utils.h" #if defined (HAVE_FFTW3) #define FFTSRC "@sc{Fftw}" #else #define FFTSRC "@sc{Fftpack}" #endif static octave_value do_fft (const octave_value_list &args, const char *fcn, int type) { octave_value retval; int nargin = args.length (); if (nargin < 1 || nargin > 3) { print_usage (); return retval; } octave_value arg = args(0); dim_vector dims = arg.dims (); octave_idx_type n_points = -1; int dim = -1; if (nargin > 1) { if (! args(1).is_empty ()) { double dval = args(1).double_value (); if (xisnan (dval)) error ("%s: NaN is invalid as the N_POINTS", fcn); else { n_points = NINTbig (dval); if (n_points < 0) error ("%s: number of points must be greater than zero", fcn); } } } if (error_state) return retval; if (nargin > 2) { double dval = args(2).double_value (); if (xisnan (dval)) error ("%s: NaN is invalid as the N_POINTS", fcn); else if (dval < 1 || dval > dims.length ()) error ("%s: invalid dimension along which to perform fft", fcn); else // to be safe, cast it back to int since dim is an int dim = NINT (dval) - 1; } if (error_state) return retval; for (octave_idx_type i = 0; i < dims.length (); i++) if (dims(i) < 0) return retval; if (dim < 0) { for (octave_idx_type i = 0; i < dims.length (); i++) if (dims(i) > 1) { dim = i; break; } // And if the first argument is scalar? if (dim < 0) dim = 1; } if (n_points < 0) n_points = dims (dim); else dims (dim) = n_points; if (dims.any_zero () || n_points == 0) return octave_value (NDArray (dims)); if (arg.is_single_type ()) { if (arg.is_real_type ()) { FloatNDArray nda = arg.float_array_value (); if (! error_state) { nda.resize (dims, 0.0); retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim)); } } else { FloatComplexNDArray cnda = arg.float_complex_array_value (); if (! error_state) { cnda.resize (dims, 0.0); retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim)); } } } else { if (arg.is_real_type ()) { NDArray nda = arg.array_value (); if (! error_state) { nda.resize (dims, 0.0); retval = (type != 0 ? nda.ifourier (dim) : nda.fourier (dim)); } } else if (arg.is_complex_type ()) { ComplexNDArray cnda = arg.complex_array_value (); if (! error_state) { cnda.resize (dims, 0.0); retval = (type != 0 ? cnda.ifourier (dim) : cnda.fourier (dim)); } } else { gripe_wrong_type_arg (fcn, arg); } } return retval; } /* %!error(fft()) %!assert(fft([]), []) %!assert(fft(zeros(10,0)), zeros(10,0)) %!assert(fft(zeros(0,10)), zeros(0,10)) %!assert(fft(0), 0) %!assert(fft(1), 1) %!assert(fft(1), 1) %!assert(fft(ones(2,2)), [2,2; 0,0]) %!assert(fft(eye(2,2)), [1,1; 1,-1]) */ DEFUN_DLD (fft, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} fft (@var{a}, @var{n}, @var{dim})\n\ Compute the FFT of @var{a} using subroutines from\n" FFTSRC ". The FFT is calculated along the first non-singleton dimension of the\n\ array. Thus if @var{a} is a matrix, @code{fft (@var{a})} computes the\n\ FFT for each column of @var{a}.\n\ \n\ If called with two arguments, @var{n} is expected to be an integer\n\ specifying the number of elements of @var{a} to use, or an empty\n\ matrix to specify that its value should be ignored. If @var{n} is\n\ larger than the dimension along which the FFT is calculated, then\n\ @var{a} is resized and padded with zeros. Otherwise, if @var{n} is\n\ smaller than the dimension along which the FFT is calculated, then\n\ @var{a} is truncated.\n\ \n\ If called with three arguments, @var{dim} is an integer specifying the\n\ dimension of the matrix along which the FFT is performed\n\ @seealso{ifft, fft2, fftn, fftw}\n\ @end deftypefn") { return do_fft (args, "fft", 0); } DEFUN_DLD (ifft, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} ifft (@var{a}, @var{n}, @var{dim})\n\ Compute the inverse FFT of @var{a} using subroutines from\n" FFTSRC ". The inverse FFT is calculated along the first non-singleton dimension\n\ of the array. Thus if @var{a} is a matrix, @code{fft (@var{a})} computes\n\ the inverse FFT for each column of @var{a}.\n\ \n\ If called with two arguments, @var{n} is expected to be an integer\n\ specifying the number of elements of @var{a} to use, or an empty\n\ matrix to specify that its value should be ignored. If @var{n} is\n\ larger than the dimension along which the inverse FFT is calculated, then\n\ @var{a} is resized and padded with zeros. Otherwise, if@var{n} is\n\ smaller than the dimension along which the inverse FFT is calculated,\n\ then @var{a} is truncated.\n\ \n\ If called with three arguments, @var{dim} is an integer specifying the\n\ dimension of the matrix along which the inverse FFT is performed\n\ @seealso{fft, ifft2, ifftn, fftw}\n\ @end deftypefn") { return do_fft (args, "ifft", 1); } /* %% fft-1.m %% %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) %% Comalco Research and Technology %% 02 May 2000 %!test %! N=64; %! n=4; %! t = 2*pi*(0:1:N-1)/N; %! s = cos(n*t); %! S = fft(s); %! %! answer = 0*t; %! answer(n+1) = N/2; %! answer(N-n+1) = N/2; %! %! assert(all( abs(S-answer) < 4*N*eps )); %% ifft-1.m %% %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) %% Comalco Research and Technology %% 02 May 2000 %!test %! N=64; %! n=7; %! t = 2*pi*(0:1:N-1)/N; %! s = cos(n*t); %! %! S = 0*t; %! S(n+1) = N/2; %! S(N-n+1) = N/2; %! %! assert(all( abs(ifft(S)-s) < 4*N*eps )); %% fft2-1.m %% %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) %% Comalco Research and Technology %% 02 May 2000 %!test %! M=16; %! N=8; %! %! m=5; %! n=3; %! %! x = 2*pi*(0:1:M-1)/M; %! y = 2*pi*(0:1:N-1)/N; %! sx = cos(m*x); %! sy = sin(n*y); %! s=kron(sx',sy); %! S = fft2(s); %! answer = kron(fft(sx)',fft(sy)); %! assert(all( all( abs(S-answer) < 4*M*N*eps ) )); %% ifft2-1.m %% %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) %% Comalco Research and Technology %% 02 May 2000 %!test %! M=12; %! N=7; %! %! m=3; %! n=2; %! %! x = 2*pi*(0:1:M-1)/M; %! y = 2*pi*(0:1:N-1)/N; %! %! sx = cos(m*x); %! sy = cos(n*y); %! %! S = kron(fft(sx)',fft(sy)); %! answer=kron(sx',sy); %! s = ifft2(S); %! %! assert(all( all( abs(s-answer) < 30*eps ) )); */ /* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: *** */