[Octave-bug-tracker] [bug #43742] ifft, ifft2, ifftn should accept 'nonsymmetric' and 'symmetric' options

[Hg200]

Follow-up Comment #5, bug #43742 (project octave):

I think this option is very important and should definitely be implemented in
Octave. Lot's of Matlab codes uses this argument and you can't easily work
around. My answer to the questions from comment #4:

1.) Does Octave test whether the input is conjugate symmetric or always use
the same FFTW routine?

IMO, Octave always uses the same complex to complex routine for the inverse
transform. If I follow the code, the calling sequence looks like this:


ifft -> do_fft at fft.cc -> ComplexNDArray::ifourier at CNDArray.cc ->
octave::fftw::ifft at oct-fftw.cc --> fftw_execute_dft at FFTW


2.) Does implementing 'symmetric' correspond to calling one of these FFTW
routines instead? 

I don't think so. It must be a transformation from complex to real.
"Symmetric" is a property of the discrete Fourier transform that says that if
x is real-valued, its discrete Fourier transform X:=DFT(x) satisfies the
property conj(X[N-n])=X[n]. For example, the following code always returns
zero:


x = rand (1, 8)
X = fft (x)
Y = conj (X([1,end:-1:2]))
X - Y


Some useful information can be found in the FFTW documentation in chapter 2.3:
"In many practical applications, the input data in[i] are purely real numbers,
in which case the DFT output satisfies the Hermitian redundancy ... It is
possible to take advantage of these circumstances in order to achieve roughly
a factor of two improvement in both speed and memory usage ... The only
differences are that the input (or output) is of type double and there are new
routines to create the plan."

Hence my conclusion:

1.) The "easiest" workaround would be to simply implement the option as: 


real(ifft(X))


2.) A much bigger patch would be to use an associated FFTW routine. I don't
know if this can be easily implemented into the current Octave code. However
the following FFTW routine seems to give the desired result for me:


fftw_complex in[NUM_POINTS];
double out[NUM_POINTS];

fftw_plan plan = fftw_plan_dft_c2r_1d (NUM_POINTS, in, out, FFTW_BACKWARD);
fftw_execute(plan);


Comments on this? Does anyone know why this was not implemented from the
beginning?

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?43742>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/