bugGNU Octave - Bugs: bug #49996, [octave forge] (signal) cheby2...

 
 

bug #49996: [octave forge] (signal) cheby2 produces zeros that are not quite complex conjugate pairs

Submitter:  None
Submitted:  Fri 06 Jan 2017 02:39:50 PM UTC
   
 
Category:  Octave Package Severity:  2 - Minor
Priority:  1 - Later Item Group:  Inaccurate Result
Status:  Confirmed Assigned to:  None
Originator Name:  Hale812 Originator Email:  -email is unavailable-
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

Fri 01 Mar 2019 06:58:35 AM UTC, comment #19: 

Oops,  Here it is.

Bill Lash <welash>
Fri 01 Mar 2019 05:46:41 AM UTC, comment #18: 

Thanks for looking at this again. There is no newcheby2.patch attached, can you try again?

Mike Miller <mtmiller>
Group Member
Fri 01 Mar 2019 04:34:10 AM UTC, comment #17: 

I just noticed that this got updated, and so I took another look at this now that cplxpair.m has been updated.  I am attaching a patch that solves, at least, the reported problem.  The change is to cheby2.m and uses cplxpair on the zero and pole vectors before calling sf2trans.  I plotted the frequency response of the resulting filter, and it looks like it meets the specification.

I'm attaching a patch called newcheby2.patch

Bill Lash <welash>
Wed 04 Apr 2018 01:58:32 AM UTC, comment #16: 

I think adding a call to cplxpair in cheby2 is the correct solution to the cheby2 issue.  I think the change to sftrans was only to address the issue with butter.  I am surprised that the butter example works now without a change to force the complex pairs during the warping process.  But if it works, that is great.

Bill Lash <welash>
Tue 03 Apr 2018 10:17:41 PM UTC, comment #15: 

Since cplxpair has been fixed in Octave, the 'butter' example mentioned in comment #4 is no longer a problem.

To sum up, the problem appears limited to the vector of zeros returned by cheby2.

In (not-yet-released) Octave 4.4, the cplxpair function forces the results that are near-conjugates, within the given (or default) tolerance, to be exact.

So if Octave 4.4 is active, we could use cplxpair within the cheby2 function to force the zeros to be exact conjugates, with some suitable tolerance.

Or some other change to ensure that the complex zeros are conjugate pairs.

Is there any reason to make a change in sftrans.m as suggested in comment #14? That is a much broader change to address a narrow problem in cheby2.

Mike Miller <mtmiller>
Group Member
Wed 11 Jan 2017 04:48:28 AM UTC, comment #14: 

I agree that the issue in butter is not with cplxpair, although I think there are some issues with cplxpair.  I was able to get the butter example to work by playing with sftrans.  In sftrans when the poles and zeros are adapted from low pass with unit cutoff to the proper type (bandpass in this case) a new set of poles and zeros are calculated. By basically forcing these to be pairs, the result will pass cplxpair.

I wrote a function that is like cplxpair, but averages the real and imaginary parts of the pairs, and called that in the bandpass case of sftrans on the b value used to calculate the poles. 

here is the patch to sftrans:

--- ../octave/signal-1.3.2/sftrans.m        2017-01-05 22:10:48.690453045 -0600
+++ sftrans.m        2017-01-07 22:46:39.951910271 -0600
@@ -141,6 +141,7 @@
 ## ----------------  -------------------------  ------------------------
       Sg = Sg * (C/(Fh-Fl))^(z-p);
       b = Sp*((Fh-Fl)/(2*C));
+b=cplxpair_eqlz(b);
       Sp = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];
       if isempty(Sz)
         Sz = zeros(1,p);


and the cplxpair_eqlz.m is


function ze = cplxpair_eqlz (z, thresh = 100*eps)

  ## interesting for testing: if nargin<2, thresh=1E-3; endif

  if isempty(z)
    ze=[];
  else
    ze = cplxpair(z); # sort complex pairs, real roots at end
    nz = length(z);
    nzrsec = 0;
    i=nz;
    while i && abs(imag(ze(i)))<thresh # determine no. of real values
      ze(i) = real(ze(i));
      nzrsec = nzrsec+1;
      i=i-1;
    endwhile
    nzsect2 = nz-nzrsec;
    if mod(nzsect2,2)~=0
      error('cplxreal: Odd number of complex values!');
    endif
    for i = 1:2:nzsect2-1
      ze(i) = (ze(i) + conj(ze(i+1)))/2;
      ze(i+1) = conj(ze(i));
    endfor
  endif

endfunction


realistically, any time a filter that produces real valued output from real valued input should have all the poles and zeros that are either real, or complex conjugate pairs.  I'm not sure if averaging the real and imaginary parts is better or worse than choosing one of the pair, and conjugating it to get the other member of the pair.

Bill Lash <welash>
Tue 10 Jan 2017 10:31:04 PM UTC, comment #13: 

Yes, the 1e-3 approximation I posted in comment #10 has the same problem in Octave as the original


>> [z, p, k] = cheby2 (6, 40, [.001, .983]);
>> cplxpair (z)
error: cplxpair: could not pair all complex numbers
error: called from
    cplxpair at line 118 column 9


But the values that SciPy computes for the same filter design call behave when pasted into Octave:


>> z2 = [(-0.99866962186808894-0.051565360057387599j),
 (-0.99928680643819134-0.037760806116408051j),
 (-0.99990441719062173-0.013825935146777633j),
 (-0.99990441719062173+0.013825935146777633j),
 (-0.99928680643819134+0.037760806116408051j),
 (-0.99866962186808894+0.051565360057387599j),
 (0.99999539574401186+0.0030345495179681286j),
 (0.99999753249436207+0.0022214871566177065j),
 (0.99999966940456275+0.0008131363754524236j),
 (0.99999966940456275-0.0008131363754524236j),
 (0.99999753249436207-0.0022214871566177065j),
 (0.99999539574401186-0.0030345495179681286j)];
>> cplxpair (z2)
ans =

  -0.999904417190622 - 0.013825935146778i
  -0.999904417190622 + 0.013825935146778i
  -0.999286806438191 - 0.037760806116408i
  -0.999286806438191 + 0.037760806116408i
  -0.998669621868089 - 0.051565360057388i
  -0.998669621868089 + 0.051565360057388i
   0.999995395744012 - 0.003034549517968i
   0.999995395744012 + 0.003034549517968i
   0.999997532494362 - 0.002221487156618i
   0.999997532494362 + 0.002221487156618i
   0.999999669404563 - 0.000813136375452i
   0.999999669404563 + 0.000813136375452i


Mike Miller <mtmiller>
Group Member
Tue 10 Jan 2017 09:46:40 PM UTC, comment #12: 

hmmm... verbatim block render fail on the website.

[z,p,k] = cheby2(6,40,[.529150859255659, 491.5561729754713269]*2/1000);

Nicholas Jankowski <nrjank>
Group Member
Tue 10 Jan 2017 09:44:31 PM UTC, comment #11: 

is the scipy version well behaved with the band edges set as in comment #1? (in octave I didn't hit the error until using quite a bit of the specified precision:

[z,p,k] = cheby2(6,40,[.529150859255659, 491.5561729754713269]*2/1000)

dunno if there's a 'nice round number' set of bug-tripping values we can use as with the butter example.

Nicholas Jankowski <nrjank>
Group Member
Tue 10 Jan 2017 08:28:21 PM UTC, comment #10: 

For those looking to fix this and with no access to Matlab (like me), or if you want a working implementation to derive code from, SciPy has free implementations of the same filter design functions, so both the outputs and the algorithms themselves can be analyzed and compared.


>>> from pprint import pprint
>>> from scipy.signal import cheby2
>>> a = cheby2 (6, 40, [.001, .983], btype='bandpass', output='zpk')
>>> (z, p, k) = (list(a[0]), list(a[1]), a[2])
>>>
>>> pprint (z)
[(-0.99866962186808894-0.051565360057387599j),
 (-0.99928680643819134-0.037760806116408051j),
 (-0.99990441719062173-0.013825935146777633j),
 (-0.99990441719062173+0.013825935146777633j),
 (-0.99928680643819134+0.037760806116408051j),
 (-0.99866962186808894+0.051565360057387599j),
 (0.99999539574401186+0.0030345495179681286j),
 (0.99999753249436207+0.0022214871566177065j),
 (0.99999966940456275+0.0008131363754524236j),
 (0.99999966940456275-0.0008131363754524236j),
 (0.99999753249436207-0.0022214871566177065j),
 (0.99999539574401186-0.0030345495179681286j)]
>>>
>>> pprint (p)
[(0.99917623666441646+0.0042929489610568049j),
 (0.99777097697199735+0.0031385196512288016j),
 (0.99696202281637791+0.0011479500050488025j),
 (0.99696202281637791-0.0011479500050488025j),
 (0.99777097697199735-0.0031385196512288016j),
 (0.99917623666441646-0.0042929489610568049j),
 (-0.98362861112405386-0.071960682314650762j),
 (-0.96149154708651674-0.051482086995079612j),
 (-0.94940131856521048-0.018602747339922173j),
 (-0.94940131856521048+0.018602747339922173j),
 (-0.96149154708651674+0.051482086995079612j),
 (-0.98362861112405386+0.071960682314650762j)]
>>>
>>> k
0.89629009720530994


Mike Miller <mtmiller>
Group Member
Tue 10 Jan 2017 08:02:13 PM UTC, comment #9: 

regarding butter, I don't think the problem is cplxpair. the error comes in the output of butter, and matlab's cplxpair also throws an error if you pass it those values. So cplxpair is at least behaving 'compatibly', but the output from butter is not in this case.

Nicholas Jankowski <nrjank>
Group Member
Sun 08 Jan 2017 04:31:32 AM UTC, comment #8: 

I took a look at the butterworth example, and it appears that in that case, forcing the conjugate pairs to be true conjugate pairs before running sftrans does not solve the issue.

I think in that case, the issue is really in cplxpair, and there is a bug report about it (https://savannah.gnu.org/bugs/?47865). 

I'm trying to follow the discussion there, but it seems like it deals with a couple of different issues.  The discussion there seems to indicate that octaves implementation multiplies the tolerance by the magnitude of the entry, but when I look at the code, It doesn't seem like it does that.  (for determining if the value is real, it does seem like the tolerance is relative, but in the case of the complex values, the check is

[v,idx] = min (abs (z(i+1:p) - conj (z(i))));
if (v > tol)
  error,,,

I think maybe changing the if to

if (v/abs(z(i)) > tol

makes some sense, but it doesn't address the problem with the butter example.

Bill Lash <welash>
Sat 07 Jan 2017 04:14:10 AM UTC, comment #7: 

In the butterworth example mentioned below, it is the poles that are not complex conjugate pairs.  This could also be an issue in the cheby2 case and we just don't have a test case that finds it. It is a little more tricky in the pole case, since the number of poles can be odd

I'll try taking a look at butter to see if there is something that can be done there.

I think that in all cases, the prototype filter that is generated will have complex conjugate poles and zeros with the differences small enough to pair, but then when the frequency is warped to generate a filter with a different cutoff or change from low-pass to high-pass, band-pass or band-stop, the small difference is amplified.

Bill Lash <welash>
Sat 07 Jan 2017 03:58:27 AM UTC, comment #6: 

The filter needs to have all the zeros in complex conjugate pairs, which is why zp2sos was complaining.  The zeros are calculated by getting equally spaced values of theta in the range of 0 to pi (not inclusive), spaced apart by pi/n (where n is the order of the filter).

The zeros are then calculated by dividing 1i by the cos(theta).  In the case where n is odd, the purely imaginary zero (at theta=pi/2) is excluded.  This should produce complex conjugate pairs.

I think the issue is that with the small numerical differences between the values of theta that produce conjugate pairs, produce small numerical differences in the cos(), which in turn produce small numerical differences in the conjugate zeros.

It is possible that it would be better that does some sort of average of the values in the conjugate pairs instead of arbitrarily choosing one of the zeros and conjugating it.

Other methods I could think of to get conjugate zeros would be to ensure that the positive and negative values of cos(theta) are of the same magnitude but opposite sign, or perhaps ensuring that the smallest value of theta(1)=pi-theta(n), theta(2)=pi-theta(n-1) ...

Bill Lash <welash>
Sat 07 Jan 2017 03:17:45 AM UTC, comment #5: 

thanks Bill. One question we both had is whether his fix is 'appropriate' mathematically. 

Is it okay to just include the conjugate zeros? Or are there legitimate cases where they need to be individually determined because the patch's assumption might overlook _?_

Nicholas Jankowski <nrjank>
Group Member
Sat 07 Jan 2017 03:13:16 AM UTC, comment #4: 

the help thread first links back to an older post with the zp2sos cplxpair error, but in that case z,p,k come from butter, not cheby2. 

here's the script that fails with butter:

>> bandHz = [100 500];
>> sampleFreq = 1001;
>> nyquistFreq = 0.5 * sampleFreq;
>> w = bandHz / nyquistFreq;
>> [z,p,k] = butter(2, w);z',p',k
ans =
   1
   1
  -1
  -1
ans =
   0.571788211592795 + 0.293698922713578i
   0.571788211592761 - 0.293698922713564i
  -0.997780778679612 - 0.002216561341880i
  -0.997780778679612 + 0.002216561341880i
k =  0.637776870938282
>> [sos,g] = zp2sos(z,p,k);
error: cplxpair: could not pair all complex numbers
error: called from
    cplxpair at line 118 column 9
    cplxreal at line 52 column 9
    zp2sos at line 78 column 10


in Matlab 2016b, the output produced is:


[z,p,k] = butter(2, w)
z =
     1
     1
    -1
    -1
p =
  0.571788211592765 + 0.293698922713572i
  0.571788211592765 - 0.293698922713572i
 -0.997780778679611 + 0.002216561341880i
 -0.997780778679611 - 0.002216561341880i
k =
   0.637776870938272
>> [sos,g] = zp2sos(z,p,k)
sos =
  Columns 1 through 4
   1.000000000000000  -2.000000000000000   1.000000000000000   1.000000000000000
   1.000000000000000   2.000000000000000   1.000000000000000   1.000000000000000
  Columns 5 through 6
  -1.143576423185530   0.413200816119565
   1.995561557359223   0.995571395446674
g =
   0.637776870938272


the patch on the thread adjusted the way zeros are calculated within cheby2, and the code ran without error except for an ordering issue. Maybe the bug needs to be broadened to butter, and maybe cheby1 and ellip as well?  is there a similarly easy test for those?

Nicholas Jankowski <nrjank>
Group Member
Sat 07 Jan 2017 03:00:20 AM UTC, comment #3: 

Here is the patch that I posted in the original discussion thread.  I have added a few comments to explain what is being done.

Basically, the change calculates only the first half of the zeros, and then conjugates those zeros and appends them.  This insures that all the zeros are complex conjugate pairs. 



(file #39392)

Bill Lash <welash>
Sat 07 Jan 2017 02:34:50 AM UTC, comment #2: 

just for reference, here's a link to the help thread with some of the technical discussion.
http://octave.1599824.n4.nabble.com/zp2sos-error-cplxpair-could-not-pair-all-complex-numbers-tp4681260.html

I've added the OP and patch author to the bug's cc line.

Nicholas Jankowski <nrjank>
Group Member
Fri 06 Jan 2017 03:48:08 PM UTC, comment #1: 

Thanks for reporting, updating summary and marking as confirmed based on mailing list discussion. Someone on the mailing list had a possible patch, hopefully that will get posted here and can be tested and merged into the signal package reasonably quickly.

Mike Miller <mtmiller>
Group Member
Fri 06 Jan 2017 02:39:50 PM UTC, original submission:  

The bug is serious preventing operation of scripts for automated filter synthesis. Sometimes, optimization code can pass speciffic values, close to Nyquist frequency, to the cheby2 function. The result is "cplxpair: could not pair all complex numbers" message and following break of the program.

Either "frequency close to Nyquist" bug should be corrrected completely, or 100% safe frequecy constraints should be specified in manuals.

fpL=0.52915085925565918945068233369966;
fpH=491.55617297547132693580351769924;
[z,p,k]=cheby2(6,40,[fpL.*2 ./1000,fpH.*2 ./1000]);
[sos, g] = zp2sos (z, p, k)


Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #46380:  newcheby2.patch added by welash (488B - text/x-patch - patch that ensures complex conjugate pairs for poles and zeros before calling sftrans())
file #39392:  cheby2.patch added by welash (870B - text/x-patch - patch to ensure complex conjugate zeros)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Updated the item)
  • -email is unavailable- added by welash (Updated the item)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by nrjank
  • -email is unavailable- added by nrjank
  •  

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

    Date Changed by Updated Field Previous Value => Replaced by
    2019-03-01 welash Attached File- Added newcheby2.patch, #46380
    2019-02-27 mtmiller Priority5 - Normal 1 - Later
        Item GroupIncorrect Result Inaccurate Result
        Operating SystemMicrosoft Windows Any
    2019-02-27 mtmiller Carbon-CopyRemoved 80942 -
    2018-05-01 mtmiller Release4.2.0 dev
    2017-08-12 jwe Summarysignal package: cheby2 produces zeros that are not quite complex conjugate pairs [octave forge] (signal) cheby2 produces zeros that are not quite complex conjugate pairs
    2017-01-07 welash Attached File- Added cheby2.patch, #39392
    2017-01-07 nrjank Carbon-Copy- Added -email is unavailable-
        Carbon-Copy- Added -email is unavailable-
    2017-01-06 mtmiller Severity3 - Normal 2 - Minor
        Item GroupSegfault, Bus Error, etc. Incorrect Result
        StatusNone Confirmed
        Summaryzp2sos error: cplxpair: could not pair all complex numbers signal package: cheby2 produces zeros that are not quite complex conjugate pairs

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code