bugGNU Octave - Bugs: bug #43742, ifft, ifft2, ifftn should accept...

 
 

bug #43742: ifft, ifft2, ifftn should accept 'nonsymmetric' and 'symmetric' options

Submitter:  Mark <mb1234>
Submitted:  Tue 02 Dec 2014 08:31:04 PM UTC
   
 
Category:  Octave Function Severity:  1 - Wish
Priority:  3 - Low Item Group:  Feature Request
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 24 Mar 2024 07:29:35 PM UTC, comment #20: 

As promised, I am attaching two patches that should cover any outstanding issues to hopefully close a 10 year old bug report. Let me briefly summarize what you get.

1.) The first patch adds two "FFTW-C2R" planners, one with single and the other with double precision. Unlike "C2C" and "R2C", "C2R" erases its input data array for rank > 1 cases. Furthermore, "C2R" only uses data from the lower half of the last dimension (in row major order) plus one element. There are auxiliary functions to handle these issues. Interface functions are necessary for the combinations of ifft, ifft2 and ifftn & single and double precision types & complex and real valued input (CNDArray, dNDArray, fCNDArray, fNDArray).

2.) The obvious prerequisite for doing an inverse Fourier transform of type "C2R" is that the input data must be a conjugate-symmetric spectrum. What happens if it is not? Well, in the rank == 1 case, the "FFTW-C2R" transformer cuts away half of the input data so that the data is treated as if it were conjugate-symmetric based on the content stored in the non-redundant part (positive frequency locations). This means that for rank == 1, asymmetries in the input data are removed implicitely.
In the generalized Nd Fourier transform (rank > 1), the symmetry is, let's say, located in "subspace" permutations. It is a cross-complex symmetry hidden in respective hyper-edges, hyper-faces, hyper-volumes etc., so things are different and much more difficult here. To eliminate asymmetries, the choice is to cross mirror the conjugate values of the lower part in each subspace, sub-subspace, etc. into the upper part. I finally found a solution to this problem, but it requires recursion. The good news is that the heap grows exponentially with rank, while the maximum recursion depth only grows linearly. I did a stack frame investigation with gdb. As expected there is a noteworthy stack consumption, but it looks okay to me. E.g. for rank = 8 with a input size array of 8^8 the "makesymmetric" stack consumption is 2kb compared to 9kb for the complete backtrace.

3.) I've added BISTs that cover more code and also trigger various internal loops and branches as there wasn't very much before.

4.) The second patch implements the SIMD macros introduced in the first patch into the existing "R2C" and "C2C" planner code. They replace the currently hard-coded padding size with sizeof()-based calculations at compile time, taking into account the size of the data type actually used. The memory alignment can now be set by the "SIMD_ALIGN_BYTES" definition and is automatically overtaken into the complete code. Note that unlike "R2C" and "C2C", "C2R" always requires an intermediate buffer, which is SIMD-aligned by default. Nevertheless, I have kept the alignment check for the input array in the "C2R" planner to maintain compatibility with the existing "R2C" and "C2C" code.

5.) Finally, a small note on the existing code. If you look at fftw::ifft(), fftw::fft(), the statement "dist = (dist < 0 ? npts : dist)" is unnecessary. However since it does not cause damage i have overtaken this into my code to be compatible with the existing code.

Questions and comments are welcome.

(file #55884, file #55885)

Hg200 <hg200>
Tue 05 Mar 2024 07:46:04 PM UTC, comment #19: 

noting that bug #65414 which was indicated as a dependency for this issue was pushed to default (so Octave v10)

Nicholas Jankowski <nrjank>
Group Member
Tue 13 Feb 2024 07:00:23 AM UTC, comment #18: 

The conclusion in comment #13 was to not use the "real part" trick to solve this problem.

I am working on a complete new patch containing "a full C2R implementation" and i am in the testing phase now.  I can release this in the next weeks. This new patch will clarify many open points:

1.) A new FFTW "C2R" planner for complex to real transformations for double and float precision. This means that "symmetric" transformations are handled natively by the correct FFTW code.
2) The FFTW code "C2R" works differently than "R2C" or "C2C". "C2R" destroys its input data. This is taken into account.
3) Asymmetries in the input data are removed by a "make symmetric" cleanup function --> Here the "real part trick" has a very big problem!!
4) New macros for SIMD alignment for better readability.
5) BISTs

Hg200 <hg200>
Tue 13 Feb 2024 04:04:24 AM UTC, comment #17: 


comment #15:

> In case anyone else would like to take a look, I already rebased the patch on a current head.
> Attaching the updated patch here.
>
>
> (file #55103)


The patch is still applies cleanly to 9.0.91.

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 13 Feb 2024 03:13:04 AM UTC, comment #16: 

Any chance to get it into 9.1 or at least 10?

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Fri 01 Sep 2023 05:00:42 PM UTC, comment #15: 

In case anyone else would like to take a look, I already rebased the patch on a current head.
Attaching the updated patch here.


(file #55103)

Markus Mützel <mmuetzel>
Group administrator
Fri 01 Sep 2023 03:46:55 PM UTC, comment #14: 

IIUC, that is point 2) in @hg200's comment #7.

Markus Mützel <mmuetzel>
Group administrator
Fri 01 Sep 2023 03:43:35 PM UTC, comment #13: 

I'm not sure if it makes sense to implement this flag in the way proposes here.
Instead, it would probably be better to setup FFTW to use a planner from the function `fftw_plan_many_dft_c2r` [1] or `fftw_plan_many_r2r` [2].

That should be done in liboctave probably around here:
https://hg.savannah.gnu.org/hgweb/octave/file/774c54278e97/liboctave/numeric/oct-fftw.cc#l247


[1]: https://www.fftw.org/fftw3_doc/Advanced-Real_002ddata-DFTs.html
[2]: https://www.fftw.org/fftw3_doc/Advanced-Real_002dto_002dreal-Transforms.html

Markus Mützel <mmuetzel>
Group administrator
Thu 31 Aug 2023 03:41:35 PM UTC, comment #12: 

What else would be needed to push the patch from comment #11 into the development version of Octave?

Guillaume <gyom>
Fri 07 Jan 2022 12:45:47 PM UTC, comment #11: 

attached is the updated patch as described in comment #10

(file #52620)

Hg200 <hg200>
Fri 31 Dec 2021 08:16:15 PM UTC, comment #10: 

I have tested fix 1.) for a few days now. There are two issues I would like to update: First, Matlab gives an error when the "symflag" is put into the "fft()" instead into the "ifft()" function. In this case, the Octave patch does not give an error message. Second, there is a "corner case" when a one-point ifft() is calculated:


ifft(1+1i,'symmetric')


In fft.cc this is handled separately, without calling FFTW. In this case, however, the real part must also be returned. It's probably also a good idea to add some BISTS to check the input validation.

I will do some further testing and will update the patch.

Hg200 <hg200>
Tue 28 Dec 2021 01:22:08 PM UTC, comment #9: 

Thank you very much for this. I agree with your recommendation: it would be great if "patch_ifft_symmetric.diff" could be applied as a first step.

Guillaume <gyom>
Sat 25 Dec 2021 03:46:58 PM UTC, comment #8: 

attaching files

(file #52565, file #52566)

Hg200 <hg200>
Sat 25 Dec 2021 03:44:45 PM UTC, comment #7: 

1.) The attached patch "patch_ifft_symmetric.diff" implements the "real(ifftx(..., 'symmetric'))" solution described in comment #6 and comment #5. The changes made are simple and limited to the "fft.cc", "fft2.cc", "fftn.cc" files. Because of the simplicity and because we don't break any architecture, my recommendation is to push the file "patch_ifft_symmetric.diff" to the mercurial repository now. If we decide to use "FFTW c2r" transformations later on, these changes are still necessary. So we don't make anything worse with this patch, but still fix the bug (see also comment #6 and comment #5).

2.) The alternative is to use the FFTW routines "fftw_plan_many_dft_c2r" and "fftw_execute_dft_c2r" according [1],[2] instead of "real(ifft())". I implemented the "c2r" routines as an example for the "ComplexNDArray" data type and attached a patch as "test_runtime_fftw_c2r_functions.diff". With this patch we can do runtime measurements and we can see where the difficulties in the implementation are: A lot of new code is added. In particular, a new method must be introduced in CNDArray.cc to transport the "symflag" information from the top software layer (fft.c) to the bottom software layer (oct-fft.cc):


NDArray
ComplexNDArray::ifouriersym (int dim) const


In oct-fft.cc there are two new methods. A new (third) FFTW planner is called to control the "c2r" transformation. This is in addition to the existing planners used for the "r2c" or "c2c" transformations.


fftw::ifftsym (const Complex *in, double *out, std::size_t npts,
               std::size_t nsamples, octave_idx_type stride,
               octave_idx_type dist)



fftw_planner::do_create_plan (const int rank, const dim_vector& dims, octave_idx_type howmany,
                              octave_idx_type stride, octave_idx_type dist,
                              const complex *in, double *out)


Regarding the runtime, we have to be careful what is compared to what. The IFFT routine is very fast and the difference is hardly noticeable when running "normal code". Are we comparing runtime when the planner is already set up, or are we including the runtime of the planner? Are we comparing runtime between processing large or small data sets? Are we comparing to "c2c" with or without taking the real part? I would say, depending on the use case, we have no advantage at worst and 30% runtime improvement at best if we compare "c2r" with "c2c" without taking the real part. If we include taking the real part, then its more.

Conclusion:
ifouriersym() needs to be also introduced in dNDArray.cc, fNDArray.cc and fCNDArray.cc. Additionally, there must be also methods for ifouriersym2() and ifouriersymn(). Where else does the code need to be changed?  Are we happy adding those new methods to the existing datatypes for just supporting a "flag"? IMO this is quite a big change.

[1]: https://www.fftw.org/fftw3_doc/Advanced-Real_002ddata-DFTs.html
[2]: http://fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html


Hg200 <hg200>
Sun 19 Dec 2021 06:06:56 PM UTC, comment #6: 

Thank you for this. I do not know when the option was introduced in Matlab but it would definitely be a good feature to have in Octave; I came across the compatibility issue again this week.

I currently use the workaround "real(ifft(X))": if the only consequence is that it is slower and uses more memory, it would still be useful to have it in place until a proper implementation is added.

Guillaume <gyom>
Sat 18 Dec 2021 05:28:23 PM UTC, comment #5: 

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?

Hg200 <hg200>
Wed 12 Feb 2020 06:27:04 PM UTC, comment #4: 

Coming across the missing 'symmetric' argument in Octave 6.0:
https://www.mathworks.com/help/matlab/ref/ifft.html#bvizm5c-1-symflag
https://www.mathworks.com/help/matlab/ref/ifft.html#bvizm5c-8

Does Octave test whether the input is conjugate symmetric or always use the same FFTW routine? Does implementing 'symmetric' correspond to calling one of these FFTW routines instead?
http://www.fftw.org/fftw3_doc/Real_002dto_002dReal-Transforms.html#Real_002dto_002dReal-Transforms

Guillaume <gyom>
Fri 16 Aug 2019 10:41:27 AM UTC, comment #3: 

Still missing in 4.4.1 (Ubuntu 19.04).

Anonymous
Sun 18 Dec 2016 08:05:34 PM UTC, comment #2: 

This feature is still missing in Octave 4.2.0.

Hartmut <hardy>
Wed 03 Dec 2014 01:24:31 AM UTC, comment #1: 

Thanks for your bug report. This change will also need to be applied to the ifft2 and ifftn functions.

Mike Miller <mtmiller>
Group Member
Tue 02 Dec 2014 08:31:04 PM UTC, original submission:  

I just came across some MATLAB code with 'symmetric' used as an argument to ifft. The documentation page describles it:

http://www.mathworks.com/help/matlab/ref/ifft.html

Current version of Octave 3.8.2 doesn't accept it:

octave:78> ifft([1 1 1 1],'symmetric')
error: invalid conversion from string to real scalar

Mark <mb1234>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #55103:  bug43742-ifft-symmetric-v3.patch added by mmuetzel (15KiB - application/octet-stream)
file #52565:  patch_ifft_symmetric.diff added by hg200 (13KiB - text/x-patch)

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by hg200 (Posted a comment)
  • -email is unavailable- added by gyom (Posted a comment)
  • -email is unavailable- added by gyom
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by mb1234 (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 21 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-03-24 hg200 Attached File- Added bug43742_implement_full_c2r_planner_release.diff, #55884
        Attached File- Added bug43742_simd_macros_for_r2c_and_c2c_release.diff, #55885
    2024-03-05 nrjank Dependencies- Depends on bugs #65414
    2024-02-13 nrjank Dependencies- bugs #65298 is dependent
    2023-09-01 mmuetzel Attached File- Added bug43742-ifft-symmetric-v3.patch, #55103
    2022-01-07 hg200 Attached File- Added patch_ifft_symmetric_update.diff, #52620
    2021-12-25 hg200 Attached File- Added patch_ifft_symmetric.diff, #52565
        Attached File- Added test_runtime_fftw_c2r_functions.diff, #52566
    2020-07-24 jwe Dependencies- bugs #58819 is dependent
    2020-02-12 gyom Carbon-Copy- Added gyom
    2019-02-27 mtmiller Priority5 - Normal 3 - Low
    2018-09-20 mtmiller Carbon-CopyRemoved 80942 -
    2018-09-20 mtmiller Severity3 - Normal 1 - Wish
        Item GroupMatlab Compatibility Feature Request
        SummaryMATLAB compatibility - ifft symmetric keyword ifft, ifft2, ifftn should accept 'nonsymmetric' and 'symmetric' options
    2018-09-20 mtmiller Dependencies- bugs #45932 is dependent
    2014-12-03 mtmiller CategoryNone Octave Function
        Item GroupNone Matlab Compatibility
        StatusNone Confirmed
        Release3.8.2 dev
        Operating SystemGNU/Linux Any

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code