bugGNU Octave - Bugs: bug #39566, Suspect interpft is wrong...

 
 

bug #39566: Suspect interpft is wrong regarding treatment of Nyquist frequency

Submitted by:  None
Submitted on:  Thu 25 Jul 2013 01:31:41 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Matlab Compatibility
Status: FixedAssigned to: Mike Miller <mtmiller>
Originator Name: Alain CochardOriginator Email: -unavailable-
Open/Closed: ClosedRelease: dev
Operating System: Any

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

Sun 04 Aug 2013 09:35:22 PM UTC, comment #4:

Fixed with the following change:

http://hg.savannah.gnu.org/hgweb/octave/rev/3a23cbde59d5

This fix will be part of the next major release of Octave, or you can pull the latest interpft function from the development source repository.

Mike Miller <mtmiller>
Project AdministratorIn charge of this item.
Thu 01 Aug 2013 08:10:29 PM UTC, comment #3:

Mike Miller writes:
> Update of bug #39566 (project octave):
>
> Status: Confirmed => In Progress
> Assigned to: None => mtmiller
>
> _____________________________________________________
>
> Follow-up Comment #2:
>
> I've got a patch for both problems you described, please test the
> attached if you are able to and let me know if this does what you
> expect.

Yes, this does what I expect.

> Also if you have access to Matlab, which it seems you do, it would
> be helpful to post some examples of interpft on complex input,
> since the Nyquist bin will not be purely real and the FFT is not
> conjugate symmetric. For example, what do you get with
>
>
> > fft (interpft ((1+j)*[1:4]', 6))

With matlab, I get:

15.0000 +15.0000i
-6.0000 - 0.0000i
-1.5000 - 1.5000i
0
-1.5000 - 1.5000i
-0.0000 - 6.0000i

Alain Cochard <cochard>
Thu 01 Aug 2013 05:27:29 AM UTC, comment #2:

I've got a patch for both problems you described, please test the attached if you are able to and let me know if this does what you expect.

I'd like to add some test cases to this change as well but haven't done those yet.

Also if you have access to Matlab, which it seems you do, it would be helpful to post some examples of interpft on complex input, since the Nyquist bin will not be purely real and the FFT is not conjugate symmetric. For example, what do you get with

(file #28731)

Mike Miller <mtmiller>
Project AdministratorIn charge of this item.
Thu 25 Jul 2013 02:25:20 PM UTC, comment #1:

Confirmed, this also affects the current development version.

Mike Miller <mtmiller>
Project AdministratorIn charge of this item.
Thu 25 Jul 2013 01:31:41 PM UTC, original submission:

As far as I understand, that is how interpft is working when the input
is a real array (I will speak only of a 1D array):

- take the discrete Fourier transform (DFT) of the array,

- zero-pad the transformed array, the zeroes being in the middle, so
as to preserve the usual symmetry properties of the DFT,

- inverse DFT of the zero-padded array.

When the length, N, of the initial vector is even, there is only one
element of the DFTed array which corresponds to the Nyquist frequency
(both the negative and positive Nyquist values are real and equal).
When zero-padding, Octave chooses to keep its associated value in the
negative range. To illustrate, simply take the DFT of the result, to
see the zero-padded array that Octave builds:

fft( interpft([1 2 3 4], 6) )'

15.00000 - 0.00000i
-3.00000 - 3.00000i
-0.00000 - 0.00000i
0.00000 - 0.00000i
-3.00000 - 0.00000i <<<< Nyquist
-3.00000 + 3.00000i

Matlab chooses instead to take half of the value associated to the
Nyquist frequency and to put it both in the positive and the negative
range:

15.0000
-3.0000 - 3.0000i
-1.5000 <<<< positive Nyquist
0
-1.5000 <<<< negative Nyquist
-3.0000 + 3.0000i

I find Matlab choice more elegant since the zero-padded array has the
usual symmetry properties: zero frequency alone, negative frequencies
complex conjugate of positive frequencies. However, the real parts of
the results of interpft are identical for Octave and Matlab, so I
would not really qualify this as a bug.

When the length, N, of the initial vector is odd, then there is a
positive and a negative Nyquist (complex conjugate of one another),
and there should be no choice for the zero-padded array. Matlab
operates in the intuitive way:

Octave, instead, puts the positive Nyquist in the negative range, next
to the negative Nyquist:

This choice is more crucial, because, in addition of the loosing of
the above mentioned symmetry, the result is different:

Matlab:

Octave:

If this is a conscious choice of Octave, it seems to me that it should
be mentioned in the documentation.

Anonymous

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #28731:  interpft.diff added by mtmiller (825B - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by cochard (Posted a comment)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by None (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 9 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Sun 04 Aug 2013 09:35:22 PM UTCmtmillerStatusIn Progress=>Fixed
      Open/ClosedOpen=>Closed
    Thu 01 Aug 2013 05:27:29 AM UTCmtmillerAttached File-=>Added interpft.diff, #28731
      StatusConfirmed=>In Progress
      Assigned toNone=>mtmiller
    Thu 25 Jul 2013 02:25:20 PM UTCmtmillerItem GroupIncorrect Result=>Matlab Compatibility
      StatusNone=>Confirmed
      Release3.6.4=>dev
      Operating SystemGNU/Linux=>Any

    Back to the top


    Powered by Savane 3.1-cleanup1