fftr now working, but very slowly

[email protected] Wed, 6 Aug 2003 15:38:47 -0300

Newsgroups	gmane.comp.gnu.octave.sources
Message-ID	<[email protected]>

Hi.

I made it work now, but it's very slow.

I want  to  modify  octave  to  make  it  use the fftw3 library, and I'm
starting with separate functions first. I expected that this real-signal
function  would work faster than octave's fft function, but of course, I
was  wrong. It's slow yet. Could it be because it's all in a dynamically
linked library, and not yet inside octave?

I never worked modifying a program's code and submitting patches, are
there any rules I should follow to contribute to octave?

bye

On Wed, Aug 06, 2003 at 11:49:13AM -0300, [email protected] wrote:
> Here is the source code of the funtion.  I'm compiling it with
> 

-- 
Nicolau Werneck <[email protected]>         9F99 25AB E47E 8724 2F71
http://cefala.org/~nwerneck                   EA40 DC23 42CE 6B76 B07F
"The secret of genius is to carry the spirit of childhood into maturity"
-- Thomas Henry Huxley

fftr.cc (text/x-c++src, 2.3 KB)

#include "/home/nwerneck/include/octave-2.1.50/octave/oct.h"
#include "/home/nwerneck/include/octave-2.1.50/octave/oct-cmplx.h"
#include "/home/nwerneck/include/octave-2.1.50/octave/lo-mappers.h"

#include "/home/nwerneck/include/fftw3.h"


//typedef size_t = unsigned int;

DEFUN_DLD(fftr, args,, "-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} fftr (@var{signal}[, @var{N}])\n\
test function for FFT of real signals using fftw3 library.
@end deftypefn")
{
// begin adapted from fft.cc 

	octave_value retval;

	int nargin = args.length();

	if (nargin < 1 || nargin > 2) {
		print_usage("fftr");
		return retval;
	}

	octave_value arg = args(0);

	int n_points = arg.rows();
	if (n_points == 1)
		n_points = arg.columns();

	if (nargin == 2) {
		double dval = args(1).double_value();
		if (xisnan(dval))
			error("fft: NaN is invalid as the N_POINTS");
		else
			n_points = NINT(dval);
	}
	if (error_state)
		return retval;

	if (n_points < 0) {
		error("fft: number of points must be greater than zero");
		return retval;
	}
	int arg_is_empty = empty_arg("fft", arg.rows(), arg.columns());

	if (arg_is_empty < 0)
		return retval;
	else if (arg_is_empty || n_points == 0)
		return octave_value(Matrix());
	if (arg.is_real_type()) {
		Matrix m = arg.matrix_value();

		if (!error_state) {
			if (m.rows() == 1)
				m.resize(1, n_points, 0.0);
			else {
				m.resize(n_points, m.columns(), 0.0);

//                      retval = m.fourier();

// begin adapted from oct-fftw.cc

				ComplexMatrix sai(m.rows() / 2 + 1, m.columns());

				double *in(m.fortran_vec());
				Complex *out(sai.fortran_vec());

				for (int i = 0; i < m.columns(); i++) {
//                                      OCTAVE_QUIT;
//                                      octave_fftw::fft(&in[npts * i],
//                                                       &out[npts * i],
//                                                       npts);

					fftw_plan ftp;
					ftp = fftw_plan_dft_r2c_1d(n_points, &in[n_points * i], reinterpret_cast < fftw_complex * >(&out[(n_points/2+1) * i]), FFTW_ESTIMATE);
					fftw_execute(ftp);
					fftw_destroy_plan(ftp);
				}

				retval = sai;
		// end adapted from oct-fftw.cc
			}

		}
	}
//  else if (arg.is_complex_type ())
//  ...
	else {
		gripe_wrong_type_arg("fft", arg);
	}
	return retval;



}