bugGNU Octave - Bugs: bug #45810, zp2sos error when calling cplxpair

 
 

bug #45810: zp2sos error when calling cplxpair

Submitter:  Charles Praplan <charprap>
Submitted:  Sun 23 Aug 2015 10:41:05 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 4.0.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 01 Sep 2015 03:24:22 PM UTC, comment #12: 
Rik <rik5>
Group administrator
Tue 01 Sep 2015 02:09:31 PM UTC, comment #11: 

Looks good to me, go ahead and push that if you're ready.

I'll repurpose this bug to the fix I still need to clean up for the signal package's helper function cplxreal.

Mike Miller <mtmiller>
Group Member
Tue 01 Sep 2015 01:02:20 PM UTC, comment #10: 

@Mike: good point.  The tolerance shouldn't be absolute--what we have now--nor should it depend only on the largest value in Z.

See the next attached changeset cpair.cset.  The tolerance is now calculated on a per pair basis (Z(i)).  While I was fixing this I decided to overhaul the function so there are a lot of other worthwhile changes in the cset. 


(file #34787)

Rik <rik5>
Group administrator
Tue 01 Sep 2015 01:26:04 AM UTC, comment #9: 

I have to make a change to the signal package as well to accommodate the original bug report, but this fix to cplxpair helps with that.

Rik, do you think this change makes the pairing too lenient? The way I read the Matlab function help, it makes it sound like the tolerance is done element-wise rather than based on the max(abs(z)).

The following example shows how this can go wrong:


octave:1> z = [2e6 + j; 2e6 - j; 1e-9 * (1 + j); 1e-9 * (1 - 2j)];
octave:2> cplxpair (z)
ans =

   1.0000e-09 - 2.0000e-09i
   1.0000e-09 + 1.0000e-09i
   2.0000e+06 - 1.0000e+00i
   2.0000e+06 + 1.0000e+00i


If the tolerance were taken from the magnitude of each element, z(3) and z(4) wouldn't be a valid complex conjugate pair.

Mike Miller <mtmiller>
Group Member
Mon 31 Aug 2015 06:52:05 PM UTC, comment #8: 

I tried it with my octave 4.0.0 and it fixes the reported problem for me.
Doug Stewart

Doug Stewart <dastew>
Mon 31 Aug 2015 04:11:25 PM UTC, comment #7: 

Yes, sorry for losing track, I will look at this later today.

Mike Miller <mtmiller>
Group Member
Mon 31 Aug 2015 03:07:40 PM UTC, comment #6: 

@Mike: Can you review my patch?  If it works I'll go ahead and commit it.

Rik <rik5>
Group administrator
Thu 27 Aug 2015 03:58:05 AM UTC, comment #5: 

@Mike: I think the problem is that the tolerance for pairing did not scale with the inputs.  If you are using numbers of order 2000, then the tolerance needs to be of order eps(2000).  Instead the function was using a fixed tolerance of 100 x eps(1).  Try the attached diff.  It works for me for the example in comment #3.  It may also fix the problem with zp2sos.

(file #34750)

Rik <rik5>
Group administrator
Sun 23 Aug 2015 08:34:10 PM UTC, comment #4: 

ML r2015b prerelease (empty lines removed):

>> z1 = [2000 + 4j; 2000 - 4j]
z1 =
   1.0e+03 *
   2.0000 + 0.0040i
   2.0000 - 0.0040i
>> cplxpair (z1)
ans =
   1.0e+03 *
   2.0000 - 0.0040i
   2.0000 + 0.0040i
>> z2 = [2000 * (1+eps) + 4j; 2000 * (1-eps) - 4j];
>> cplxpair (z2)
ans =
   1.0e+03 *
   2.0000 - 0.0040i
   2.0000 + 0.0040i


Philip Nienhuis <philipnienhuis>
Group Member
Sun 23 Aug 2015 07:19:42 PM UTC, comment #3: 

I think there is a compatibility bug with cplxpair here, but maybe also with zp2sos once that is fixed.

Example:

>>> z1 = [2000 + 4j; 2000 - 4j];
>>> cplxpair (z1)
ans =

   2000 -    4i
   2000 +    4i

>>> z2 = [2000 * (1+eps) + 4j; 2000 * (1-eps) - 4j];
>>> cplxpair (z2)
error: cplxpair: could not pair all complex numbers
error: called from
    cplxpair at line 128 column 9


Can OP or someone try this in Matlab to verify that it should work?

My reading of Matlab's help on cplxpair says that this should work because the tolerance is taken relative to the magnitude of each value.

Mike Miller <mtmiller>
Group Member
Sun 23 Aug 2015 07:09:11 PM UTC, comment #2: 

Is there a simple example that can be used to demonstrate the issue with cplxpair?  The first step is to determine whether this is a problem with code in the signal package or in core Octave.

Rik <rik5>
Group administrator
Sun 23 Aug 2015 11:52:31 AM UTC, comment #1: 

this work around will help you until we get it fixed.

pkg load signal
N=2;
Wp=[5976, 6605];
[Z,P,K] = butter(N,Wp,'stop','s')
p1=round(1e11*P)/1e11 # one less digit of accuracy.
[SOS,G] = zp2sos(Z,p1,K)

Doug Stewart <dastew>
Sun 23 Aug 2015 10:41:05 AM UTC, original submission:  

Hi,
The following code produces an error in the function zp2sos. Note that it is working with Matlab and also if the filter type is ellip.

Thank you



% NOT WORKING

pkg load signal
N=2;
Wp=[5976, 6605];
[Z,P,K] = butter(N,Wp,'stop','s')
[SOS,G] = zp2sos(Z,P,K)


% WORKING

[Z,P,K] = ellip(N,3,40,Wp,'stop','s')
[SOS,G] = zp2sos(Z,P,K)



Charles Praplan <charprap>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #34787:  cpair.cset added by rik5 (6KiB - application/octet-stream)
file #34750:  cplxpair.diff added by rik5 (922B - text/x-diff)

 

Carbon-Copy List
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by charprap (Submitted the item)
  • -email is unavailable- added by charprap
  •  

    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 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2017-01-27 rik5 Dependencies- bugs #50124 is dependent
    2016-05-06 rik5 Dependencies- bugs #47865 is dependent
    2015-09-21 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2015-09-01 rik5 Attached File- Added cpair.cset, #34787
    2015-09-01 mtmiller StatusPatch Submitted In Progress
        Operating SystemMicrosoft Windows Any
    2015-08-30 rik5 StatusConfirmed Patch Submitted
    2015-08-27 rik5 Attached File- Added cplxpair.diff, #34750
        Item GroupSegfault, Bus Error, etc. Inaccurate Result
        StatusNone Confirmed
    2015-08-23 charprap Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code