bugGNU Octave - Bugs: bug #45932, ifft(fft(x)) is not real for x real

 
 

bug #45932: ifft(fft(x)) is not real for x real

Submitter:  Marco Caliari <caliari>
Submitted:  Fri 11 Sep 2015 02:42:12 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Duplicate Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 20 Sep 2018 03:39:27 PM UTC, comment #7: 

Yes, I think that is the 'c2r' FFTW planner, that is what should be implemented for bug #43742.

http://www.fftw.org/fftw3_doc/One_002dDimensional-DFTs-of-Real-Data.html

Mike Miller <mtmiller>
Group Member
Thu 20 Sep 2018 06:35:36 AM UTC, comment #6: 

@Mike: I'm pretty sure that FFTW has a special method for conjugate symmetric input. And such a method is of course faster than the normal one. The 'symmetric' option is used by Matlab to force that special method in case the input is not exactly conjugate symmetric because of rounding errors. It is not an option that simply removes the imaginary parts from the output. Even if it could be implemented as such in octave.

Marco Caliari <caliari>
Group Member
Thu 20 Sep 2018 01:15:10 AM UTC, comment #5: 

It's possible that Matlab is doing a test for whether the input is conjugate symmetric and using a different transform method, which Octave is not doing. Is it worth it or is developer time better spent on implementing the 'symmetric' option to force it? Maybe both at the same time. If someone thinks it's important to do automatic symmetry detection this can be reopened, but not a lot has happened here in 3 years.

Mike Miller <mtmiller>
Group Member
Thu 20 Sep 2018 12:59:41 AM UTC, comment #4: 

It seems to me that this bug is exactly bug #43742.

In the normal case, ifft(fft(rand(128,1))) returns a complex vector. If the user wants to force the conjugate symmetric inverse transform, it should be done with ifft(fft(rand(128,1)),'symmetric'). That is precisely bug #43742. I don't think there is any other appropriate way to do this. Am I missing something? Closing as a duplicate.

Mike Miller <mtmiller>
Group Member
Thu 25 Jan 2018 07:43:10 AM UTC, comment #3: 

I do not think that removing the imag part if small is safe. The problem in fact is the choice of K. But, there is no heuristic. In xhat of the original example, xhat(1+i) is exactly equal to conj(xhat(n+1-i)). This should be the check and the only case in which to apply the zeroing of the imaginary part.

As a second step, I think to remember that there is a proper function in FFTW to manage this case, which does not compute at all the imaginary part, but Octave does not use it (at least, when I submitted the report).

Marco Caliari <caliari>
Group Member
Wed 24 Jan 2018 06:11:12 PM UTC, comment #2: 

We can check if max(abs(imag()) is smaller then K*eps, and zero imag() if the answer is yes.
But how can we set K? Is it depend on vector/matrix size?

Avinoam Kalma <avinoam>
Group Member
Sat 12 Sep 2015 12:40:57 AM UTC, comment #1: 

No, there is no detection that the input is conjugate symmetric. If the (ifft) input array is complex, then a complex-to-complex inverse FFT is done and the output is complex.

It might make sense to add some heuristics to detect some of the symmetric special cases, if not to choose a different function then to at least clean up the round-off errors in the output.

Also note the roundoff error in double vs single:


>> max (abs (imag (ifft (fft (rand (101, 1, "double"))))))
ans =    1.2884e-16
>> max (abs (imag (ifft (fft (rand (101, 1, "single"))))))
ans =    6.0432e-08


Mike Miller <mtmiller>
Group Member
Fri 11 Sep 2015 02:42:12 PM UTC, original submission:  

If x is a real vector, say


n = 101;
x = rand(n,1);


then


xhat = fft(x);


has the property xhat(1+i)=conj(xhat(n+1-i)). In this, it seems to me that Octave uses the rfftw transform of FFTW. But if we now compute


ifft(xhat)


we get a complex vector with spurious imaginary parts. For comparison, Matlab gives a real vector. Would'n it possibile to use an (inverse) rfftw transform of FFTW? Or, alternatively, check the input and return only the real part? Notice that this report is related to #43742.

Marco

Marco Caliari <caliari>
Group Member

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Carbon-Copy List
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by caliari (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 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2024-02-13 nrjank Dependencies- bugs #65298 is dependent
    2018-09-20 mtmiller Carbon-CopyRemoved 80942 -
    2018-09-20 mtmiller StatusConfirmed Duplicate
        Open/ClosedOpen Closed
        Dependencies- Depends on bugs #43742
    2015-09-12 mtmiller StatusNone Confirmed
        Release4.0.0 dev

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code