bugGNU Octave - Bugs: bug #43673, fft() gives incorrect result for...

 
 

bug #43673: fft() gives incorrect result for 3-dim array with singleton 3rd dimension

Submitter:  Eddy <count>
Submitted:  Sun 23 Nov 2014 02:22:38 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed 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

Tue 29 Nov 2016 09:45:58 PM UTC, comment #15: 

Here is a complete patch that also fix fourier forward and inverse transforms in *NDArray.cc.

This patch could be divided into three parts:

(1) The first part is the same as file #39079, see previous comment.

(2) The second part is the fix to *NDArray.cc. Which also deal with the "DIM > ndims (X)" problem.

(3) The third part is the removal of `if (n_points == 1)` in fft.cc. Because with (1)(2), special treatment of `n_points == 1` is not necessary.

(file #39105)

Eddy <count>
Sun 27 Nov 2016 06:40:19 AM UTC, comment #14: 

This bug is not fully solved.

In fft (X, N, DIM), when DIM > ndims (X), fft should still produce reasonable result instead of casting an error.

For example:


> fft(1, 3, 3)
error: fft: DIM must be a valid dimension along which to perform FFT


Expected result:


ans =

ans(:,:,1) =  1
ans(:,:,2) =  1
ans(:,:,3) =  1


In general (when DIM > ndims (X)), the answer is repmat (X, [ones(1,DIM-1) N]). Which can be understood as 1-padding to dimensions of X to DIM, then 0 padding or truncation X to length N in dimension DIM, then do usual fft.

I have gave the use case in comment #3. And this extension conform to the behaviour of MATLAB also.

I attatched a new patch. The fix works by 1-padding dimension of X to DIM.

(file #39079)

Eddy <count>
Mon 21 Nov 2016 10:11:17 PM UTC, comment #13: 

I fixed this problem for bug #49026.  It is now part of the 4.2.0 release.  Closing report.

Rik <rik5>
Group administrator
Sat 10 Jan 2015 12:34:39 AM UTC, comment #12: 

Okay, I understand what you are doing now.  I think this is a fine thing.  You don't need to join -email is unavailable- unless you think you will become a regular contributor and want to keep abreast of what others are doing.

Some helpful instructions for contributing a changeset are at https://www.gnu.org/software/octave/doc/interpreter/Contributing-Guidelines.html.  You can look at the sub-chapters on preparing a changeset, and the ones on patching C++ code.  In general, just try and follow the existing conventions as regard to indent spacing, etc.

When you have the full patch just upload it here and I can review it and commit it.  We generally check in code with a "USER NAME <email@address>" comment.  It acknowledges those who contribute, and also gives us someone to contact if things do eventually go wrong with the patch.  But if you wish to remain anonymous I can commit it under my credentials.


Rik <rik5>
Group administrator
Fri 09 Jan 2015 08:23:44 PM UTC, comment #11: 

Thank you for your reply! I just have no experience of contributing code to a project, and don't know how the whole process flow go. (Should I join the mailing list  maintainers@octave.org?)

The dimensions resized (as you pointed out) because I want to deal with these cases (and the expected results):


> a = [1 2;3 4];
> fftlen=0; fft( a, fftlen, 3)
ans = [](2x2x0)
> fftlen=1; fft( a, fftlen, 3)
ans =

   1   2
   3   4

> fftlen=2; fft( a, fftlen, 3)
ans =

ans(:,:,1) =

   1   2
   3   4

ans(:,:,2) =

   1   2
   3   4


Note the "a" in the code could come from "a=b(2,2,1:1)", so this 2-dims to 3-dims transform is not weird.

So, except for fftlen==1, I have to resize the dimensions, since the result has dimensions like that.

Of course, these are trival cases of Fourier Transform. Here I choose a way that fix the bug with least code, instead of treat these cases specially. Not 100% sure if this is the right thing to do.
Oh, Maybe I shoud mension these concerns when providing the patch.


For dNDarray.cc, we don't need to worry the resize of dimensions because there is no the "N" parameter (fft length) there


-- Loadable Function:  fft (X, N, DIM)
v.s
ComplexNDArray
NDArray::fourier (int dim) const


So for the case of "dim >= dims.length ()", the dNDarray.cc code is always doing "fft (X, 1, DIM)", which is the special case mensioned above (fftlen==1). Return the result immediately is clearly more understandable than fix the index range logic in NDArray::fourier().


Eddy <count>
Tue 06 Jan 2015 05:26:01 PM UTC, comment #10: 

Sorry, we're a skeleton crew that takes care of Octave.  The patch generally looks good.  One question, in fft.cc if the desired dim is greater than the dimensions of the existing object it resizes the dimensions of the object and then passes it down to liboctave for the rest of the FFT processing.


     if (dims(i) < 0)
       return retval;

+  if (dim >= dims.length ())
+    {
+      dims = dims.redim (dim + 1);
+    }
+


Why not just short-circuit at this point and return the original argument?


+  if (dim >= dims.length ())
+    {
+      retval = arg;
+      return retval;
+    }


I agree that you also have to add the same code in dNDarray.cc and its relatives because people can link directly against liboctave if they want too without using the interpreter.


+  if (dim >= dv.length ())
+    {
+      ComplexNDArray retval (*this);
+      return retval;
+    }



Rik <rik5>
Group administrator
Tue 06 Jan 2015 10:21:08 AM UTC, comment #9: 

Any comment about the fix? If you feel that kind of fix is OK, I may try fix other related functions also.

Eddy <count>
Tue 16 Dec 2014 02:17:52 PM UTC, comment #8: 

I wrote a patch that fix the bug for double array when have fftw.

There are other FFT functions need to be fixed. (for different array types, with or without fftw, inverse or forward transform)


Tests that related to this bug

a=reshape(1:2*2*2, 2,2,2);
assert (fft (a, 1, 3), a(:,:,1))
assert (fft (a, 0, 3), a(:,:,1:0))
assert (fft (a(:,:,1), 1, 3), a(:,:,1))
assert (fft (a(:,:,1), 0, 3), a(:,:,1:0))
assert (fft (ones(2,2,1,2), 1, 4), ones (2,2))
assert (fft (1, 3, 3), ones (1,1,3))


(file #32681)

Eddy <count>
Sun 14 Dec 2014 08:04:14 AM UTC, comment #7: 

The function NDArray::fourier is still using uninitialized memory in the example case

> fft( zeros(2,2,2), 1, 3 )


zeros(2,2,2) is been copied(during resize) and dimension squeezed to zeros(2,2) (liboctave/array/Array.cc:1033 Array<T>::resize).

Then in liboctave/array/dNDArray.cc (and probably also exist in other xxNDArray.cc) the check in line 104
  if (dim > dv.length () || dim < 0)
fail to recognize this case (dim==2, dv.length ()==2). I guess the author want something like
  if (dim >= dv.length () || dim < 0)

Then the access (line 108)
  octave_idx_type n = dv(dim);
conceptually should be an out of bound access. In my debug, it returns n = 32 (totally wrong).

As a consequence, nloop = 0, so the main loop in line 123
  for (octave_idx_type k = 0; k < nloop; k++)
does not even run.

This finally lead to the space allocate in line 119,120
  ComplexNDArray retval (dv);
  Complex *out (retval.fortran_vec ());
been returned without initialization.


The inverse Fourier interface NDArray::ifourier returns a correct final result because its space allocation is (line 150)
  ComplexNDArray retval (*this);
So it is already been initialized (despite of nloop==0).


BTW: I see here a potential speed improvement of NDArray::ifourier, it is using in-place FFT, which is slow than the "out-place" FFT in NDArray::fourier (maybe by a factor of 2).

Eddy <count>
Sat 13 Dec 2014 06:57:42 PM UTC, comment #6: 

It looks like been fixed for the case

> fft( zeros(2,2,2), 1, 3 )


But

> fft( ones(2,2,2), 1, 3 )

ans =

   0   0
   0   0
Which is a wrong result. (expecting ones(2,2))


Also

> fft( zeros(2,2,1), 1, 3 )

still does not work as expected.

Eddy <count>
Fri 05 Dec 2014 08:11:10 PM UTC, comment #5: 

Can the reporter try a development version later than this changeset:


changeset:   19381:af41e41ad28e
user:        Kai T. Ohlhus <k.ohlhus@gmail.com>
date:        Fri Dec 05 13:08:36 2014 +0100
summary:     replace oct-mem.h inline indirections by standard function calls.


I see the problems before this cset, but after it I get the correct 2x2 matrix of zeros.

Rik <rik5>
Group administrator
Mon 24 Nov 2014 07:57:46 PM UTC, comment #4: 

The FT of length 1 is trivial, the FFT not quite so.
Anyhow, the FFT function should handle both your cases
appropriately (and some others as well).

This is a helpful report.

Michael Godfrey <godfrey>
Group Member
Mon 24 Nov 2014 02:29:34 PM UTC, comment #3: 

yeah, FFT of length 1 is trivial. But my code could be more compact without this bug...

Code in my mind:

for k = 1 : n
  wA = fft(A(:,:,1:k), n, 3);
  ...some code that involve wA
end


Currently I have to treat both the case k==1 and n==1 specially...

Eddy <count>
Mon 24 Nov 2014 12:40:55 PM UTC, comment #2: 

By the way, computing an FFT on length 1 is not
typically very useful, ut it should be andled correctly.
Other 3 argument cases are not quite right:
octave:25> fft( zeros(2,2,2), 0,3 )
ans = [](2x2x0)

So, the length 1 case needs fixing and argument checks on
the three argument case need improving, and tests for
three argument cases are needed.

Finally, last night I found a case that ended in a seg fault,
but I cannot reproduce it now. It was not written to the
history due to the seg fault.

I will report it if I can remember the arguments.

Michael Godfrey <godfrey>
Group Member
Mon 24 Nov 2014 10:54:25 AM UTC, comment #1: 

There appears to be a problem with the three argument logic.
In this example line 90 of fft.cc will set dim = 0, which
does not seem right.

Also, there are no tests for the 3 argument cases.

Michael Godfrey <godfrey>
Group Member
Sun 23 Nov 2014 02:22:38 PM UTC, original submission:  

like this:

octave:7> fft( zeros(2,2,2), 1, 3 )
ans =

    0.0000e+00 -         NaNi   1.8354e-316 +  0.0000e+00i
   2.1220e-314 + 1.6347e-316i    0.0000e+00 + 9.3214e-315i

It should return zeros(2,2).
(interestingly, ifft() gives correct result)
Bug also exist in version 3.6.2 and 3.8.1.


Also a closely related Matlab Compatibility problem:

octave:11> fft( zeros(2,2,1), 1, 3 )
error: fft: DIM must be a valid dimension along which to perform FFT

Expected result is return zeros(2,2,1) it self.
ifft() has this problem too.

Eddy <count>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #39105:  bug_43673_v3.patch added by count (15KiB - text/x-patch - fft,ifft,fourier,ifourier: Give reasonable output when DIM is greater than dimension of X.)
file #39079:  bug_43673_v2.patch added by count (2KiB - text/x-patch - fft: Give reasonable output when DIM is greater than dimension of X.)
file #32681:  fft_bug43673_pt1.diff added by count (2KiB - text/x-patch - fix the fft() bug for double Array when have fftw)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by count (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 9 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-11-29 count Attached File- Added bug_43673_v3.patch, #39105
    2016-11-27 count Attached File- Added bug_43673_v2.patch, #39079
    2016-11-21 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2016-11-21 rik5 Dependencies- Depends on bugs #49026
    2015-01-10 rik5 StatusPatch Submitted In Progress
    2014-12-29 rik5 StatusConfirmed Patch Submitted
    2014-12-16 count Attached File- Added fft_bug43673_pt1.diff, #32681
    2014-11-25 mtmiller StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code