zp2sos error: cplxpair: could not pair all complex numbers

classic Classic list List threaded Threaded
23 messages Options
12
Reply | Threaded
Open this post in threaded view
|

zp2sos error: cplxpair: could not pair all complex numbers

hale812
This post was updated on .
Again, the problem exactly as described 5 years ago. And nothing was fixed:
http://octave.1599824.n4.nabble.com/Signal-zp2sos-and-cplxpair-problem-td4597918.html

Trying to calculate 6th order Chebyshev filter, getting an error out of the blue.

The code:
conf.stages=6; %Stages Chebyshev type2 filter
conf.reject=40; %-dB out-of-band rejection
conf.fsamp=1000; %1000Hz sampling frequency
fpL=0.52915085925565918945068233369966; %Hz singular rejection frequency (pole) below the band
fpH=491.55617297547132693580351769924; %Hz singular rejection frequency (pole) above the band
%[fpL.*2 ./conf.fsamp,fpH.*2 ./conf.fsamp]=[0.00107164775883543,   0.00163819175474634]

[z,p,k]=cheby2(conf.stages,conf.reject,[fpL.*2 ./conf.fsamp,fpH.*2 ./conf.fsamp]);

error: cplxpair: could not pair all complex numbers
error: called from
    cplxpair at line 128 column 9
    cplxreal at line 52 column 9
    zp2sos at line 77 column 10

These (fpL, fpH) numbers look weird, but that is what fsolve has suggested, when optimizing the filter. Technically, they are completely legit and should not crash.


This is really an awful kind of problem. Octave is becoming popular among engineers, but such fundamental problems, like failing signal toolbox, or crashing when computing 32bit numbers in long loops, reduce it's reputation to nothing. I just can not finish my work without rewriting everything in Fortran or C.

5 years guys! The problem there is older than 5 years and nobody done a thing about it.


>> ver
GNU Octave Version: 4.0.2
Operating System: MINGW32_NT-6.1 Windows 7 Service Pack 1 i686
communications-1.2.1
control-3.0.0
signal-1.3.2
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Mike Miller-4
On Wed, Jan 04, 2017 at 21:30:14 -0800, hale812 wrote:

> Trying to calculate 6th order Chebyshev filter, getting an error out of the
> blue.
>
> The code:
> conf.stages=6; %Stages Chebyshev type2 filter
> conf.reject=40; %-dB out-of-band rejection
> conf.fsamp=1000; %1000Hz sampling frequency
> fpL=0.52915085925565918945068233369966; %Hz singular rejection frequency
> (pole) below the band
> fpH=491.55617297547132693580351769924; %Hz singular rejection frequency
> (pole) above the band
> %[fpL.*2 ./conf.fsamp,fpH.*2 ./conf.fsamp]=[0.00107164775883543,  
> 0.00163819175474634]
>
> [z,p,k]=cheby2(conf.stages,conf.reject,[fpL.*2 ./conf.fsamp,fpH.*2
> ./conf.fsamp]);
>
> error: cplxpair: could not pair all complex numbers
> error: called from
>     cplxpair at line 128 column 9
>     cplxreal at line 52 column 9
>     zp2sos at line 77 column 10

I'm guessing you omitted the last line in your script, something like

  [sos, g] = zp2sos (z, p, k);

I get the same error as you with the latest version of Octave. There's
nothing wrong with your numbers or with what you're trying to do. We
could use some help from you to make Octave better instead of a
complaint and an expectation that it somehow gets better on its own.

Here is the reason for the error

  >> [z, p, k] = cheby2 (6, 40, [fpL, fpH] .* 2 ./ 1000);
  >> [~, idx] = sort (real (z));
  >> z = z(idx).'
  z =

    -0.999905676503125 + 0.013734558486826i
    -0.999905676503125 - 0.013734558486826i
    -0.999296201278521 + 0.037511359750325i
    -0.999296201278521 - 0.037511359750325i
    -0.998687143937134 + 0.051224881988043i
    -0.998687143937134 - 0.051224881988043i
     0.999994843220576 + 0.003211468862565i
     0.999994843220576 - 0.003211468862565i
     0.999997236381364 + 0.002351006089708i
     0.999997236381364 - 0.002351006089708i
     0.999999629730722 + 0.000860545419713i
     0.999999629730722 - 0.000860545419685i

Now as you can see by the sorted pairs of zeros, the first 5 pairs are
within tolerance of each other. The last pair is not (the imaginary
values differ by about 2.8e-14), which is why the call to cplxpair
fails.

Is the error here due to cplxpair (should the default tolerance be
looser)? Or is the error due to the list of zeros returned by cheby2
(should the computed pairs be within tolerance to begin with)?

If you can help identify what Matlab does in this situation, and what
should be fixed and where, then maybe something can be done about it.

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Mike Miller-4
On Wed, Jan 04, 2017 at 22:59:29 -0800, Mike Miller wrote:

> Here is the reason for the error
>
>   >> [z, p, k] = cheby2 (6, 40, [fpL, fpH] .* 2 ./ 1000);
>   >> [~, idx] = sort (real (z));
>   >> z = z(idx).'
>   z =
>
>     -0.999905676503125 + 0.013734558486826i
>     -0.999905676503125 - 0.013734558486826i
>     -0.999296201278521 + 0.037511359750325i
>     -0.999296201278521 - 0.037511359750325i
>     -0.998687143937134 + 0.051224881988043i
>     -0.998687143937134 - 0.051224881988043i
>      0.999994843220576 + 0.003211468862565i
>      0.999994843220576 - 0.003211468862565i
>      0.999997236381364 + 0.002351006089708i
>      0.999997236381364 - 0.002351006089708i
>      0.999999629730722 + 0.000860545419713i
>      0.999999629730722 - 0.000860545419685i
>
> Now as you can see by the sorted pairs of zeros, the first 5 pairs are
> within tolerance of each other. The last pair is not (the imaginary
> values differ by about 2.8e-14), which is why the call to cplxpair
> fails.

This was simply looking into the stated problem and show that you can
use Octave's interpreter and debugging capabilities to see why the
function fails.

As a filtering problem, as was suggested in the earlier post you
referred to, you may be trying to construct a filter with too tight a
constraint. Approximating as

  cheby2 (6, 40, [.001, .983])

you are trying to design a bandpass filter that is incredibly close to
an all-pass filter, especially at the low end.

If what you want is a lowpass at 491.5 Hz, plus DC rejection, maybe
design a lowpass filter instead with

  cheby2 (6, 40, .983)

and cascade with a separate highpass filter for tight DC rejection?

Just suggestions, and all of this is not to say that there aren't still
possible bugs in any of the bilinear, cheby2, or cplxpair functions.

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

hale812
This post was updated on .
Of course, I forgot to copypaste
[sos,gain] = zp2sos(z,p,k);
But it is in the header, anyway.

So yes, the frequency is close, but not critically.

Definitely this has little physical meaning.

But it has meaning for the optimizer (fsolve) that does not care about physics. Only about numbers.

The goal was to optimize the filter stack (DSP blocks) to certain passband. AFAIK, there is no explicit formula, so I was using numerical solver of arbitrary nonlinear cost function. And the solver just happen to try some peculiar numbers from time to time. There is no crime in that, am I right?

And this bug just does not allow me to optimize higher order circuits, failing at some step with the error mentioned.

It is just a common sense, that the signal toolbox SHOULD work from 0 to Nyquist, at least!
And from my opinion, it should work in all the Euclidean space for combined computations with "negative time" and "negative frequency", which happens quite often. That depends on return value, of course.
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank
In reply to this post by Mike Miller-4
> If you can help identify what Matlab does in this situation, and what
> should be fixed and where, then maybe something can be done about it.

Running the same script in Matlab 2016b produces no errors and the
following output.
Note Octave output from cheby2 is a row vector, Matlab is a column. I
don't know if this matters. Matlab output follows. Octave 4.2.0,
signal pkg v.1.3.2 ouput shown below that.

-------------------------------------
Matlab 2016b

>> conf.stages=6;
>> conf.reject=40;
>> conf.fsamp=1000;
>> fpL=0.52915085925565918945068233369966;
>> fpH=491.55617297547132693580351769924;
>> [z,p,k]=cheby2(conf.stages,conf.reject,[fpL.*2 ./conf.fsamp,fpH.*2./conf.fsamp])
z =
 -0.999905676503125 + 0.013734558486828i
 -0.999905676503125 - 0.013734558486828i
 -0.998687143937134 + 0.051224881988044i
 -0.998687143937134 - 0.051224881988044i
 -0.999296201278521 + 0.037511359750324i
 -0.999296201278521 - 0.037511359750324i
  0.999999629730721 + 0.000860545419704i
  0.999999629730721 - 0.000860545419704i
  0.999994843220577 + 0.003211468862567i
  0.999994843220577 - 0.003211468862567i
  0.999997236381365 + 0.002351006089707i
  0.999997236381365 - 0.002351006089707i
p =
 -0.983753121609877 + 0.071492541325641i
 -0.983753121609877 - 0.071492541325641i
 -0.961750098624975 + 0.051154896580545i
 -0.961750098624975 - 0.051154896580545i
 -0.949728425507133 + 0.018486078667500i
 -0.949728425507133 - 0.018486078667500i
  0.999127670635969 + 0.004543010951549i
  0.999127670635969 - 0.004543010951549i
  0.997640877892379 + 0.003321080212941i
  0.997640877892379 - 0.003321080212941i
  0.996785136655982 + 0.001214672037223i
  0.996785136655982 - 0.001214672037223i
k =
   0.896584733013338
>> [~, idx] = sort (real (z));
>> z = z(idx)
z =
 -0.999905676503125 + 0.013734558486828i
 -0.999905676503125 - 0.013734558486828i
 -0.999296201278521 + 0.037511359750324i
 -0.999296201278521 - 0.037511359750324i
 -0.998687143937134 + 0.051224881988044i
 -0.998687143937134 - 0.051224881988044i
  0.999994843220577 + 0.003211468862567i
  0.999994843220577 - 0.003211468862567i
  0.999997236381365 + 0.002351006089707i
  0.999997236381365 - 0.002351006089707i
  0.999999629730721 + 0.000860545419704i
  0.999999629730721 - 0.000860545419704i

--------------------------------------------------------------
Octave 4.2.0

>>pkg load signal
>> conf.stages=6;
>> conf.reject=40;
>> conf.fsamp=1000;
>> fpL=0.52915085925565918945068233369966;
>> fpH=491.55617297547132693580351769924;
>> [z,p,k]=cheby2(conf.stages,conf.reject,[fpL.*2 ./conf.fsamp,fpH.*2./conf.fsamp]);
>> z'
ans =
  -0.998687143937134 - 0.051224881988043i
  -0.999296201278521 - 0.037511359750325i
  -0.999905676503125 - 0.013734558486826i
   0.999999629730722 - 0.000860545419713i
   0.999997236381364 - 0.002351006089708i
   0.999994843220576 - 0.003211468862565i
   0.999994843220576 + 0.003211468862565i
   0.999997236381364 + 0.002351006089708i
   0.999999629730722 + 0.000860545419685i
  -0.999905676503125 + 0.013734558486826i
  -0.999296201278521 + 0.037511359750325i
  -0.998687143937134 + 0.051224881988043i
>> p'
ans =
   0.999127670635973 - 0.004543010951549i
   0.997640877892371 - 0.003321080212942i
   0.996785136655987 - 0.001214672037224i
   0.996785136655987 + 0.001214672037222i
   0.997640877892375 + 0.003321080212938i
   0.999127670635973 + 0.004543010951549i
  -0.983753121609872 + 0.071492541325644i
  -0.961750098624987 + 0.051154896580549i
  -0.949728425507125 + 0.018486078667483i
  -0.949728425507125 - 0.018486078667483i
  -0.961750098624987 - 0.051154896580549i
  -0.983753121609872 - 0.071492541325644i
>> k
k =  0.896584733013339
>> [~, idx] = sort (real (z));
>> z = z(idx)'
z =
  -0.999905676503125 - 0.013734558486826i
  -0.999905676503125 + 0.013734558486826i
  -0.999296201278521 - 0.037511359750325i
  -0.999296201278521 + 0.037511359750325i
  -0.998687143937134 - 0.051224881988043i
  -0.998687143937134 + 0.051224881988043i
   0.999994843220576 - 0.003211468862565i
   0.999994843220576 + 0.003211468862565i
   0.999997236381364 - 0.002351006089708i
   0.999997236381364 + 0.002351006089708i
   0.999999629730722 - 0.000860545419713i
   0.999999629730722 + 0.000860545419685i

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank
On Thu, Jan 5, 2017 at 12:30 PM, Nicholas Jankowski <[hidden email]> wrote:
>> If you can help identify what Matlab does in this situation, and what
>> should be fixed and where, then maybe something can be done about it.
>
> Running the same script in Matlab 2016b produces no errors and the
> following output.
> Note Octave output from cheby2 is a row vector, Matlab is a column. I
> don't know if this matters. Matlab output follows. Octave 4.2.0,
> signal pkg v.1.3.2 ouput shown below that.

Oops, I left off the important part of my scripts as well.

-----------------------
Matlab 2016b:
>>  [sos, g] = zp2sos (z, p, k)
sos =
  Columns 1 through 5
   1.000000000000000   1.999811353006249   0.999999999999999
1.000000000000000   1.899456851014267
   1.000000000000000   1.998592402557042   1.000000000000000
1.000000000000000   1.923500197249950
   1.000000000000000   1.997374287874268   0.999999999999998
1.000000000000000   1.967506243219755
   1.000000000000000  -1.999999259461443   0.999999999999999
1.000000000000000  -1.993570273311965
   1.000000000000000  -1.999994472762730   1.000000000000001
1.000000000000000  -1.995281755784758
   1.000000000000000  -1.999989686441154   1.000000000000001
1.000000000000000  -1.998255341271938
  Column 6
   0.902325817320760
   0.927580075649315
   0.972881387742377
   0.993582084086443
   0.995298350815657
   0.998276741178964
g =
   0.896584733013338

-----------
Octave 4.2.0:

>> [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 77 column 10

%% just to verify it's not as simple as the orientation of z and p:

>> [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 77 column 10

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Mike Miller-4
In reply to this post by nrjank
On Thu, Jan 05, 2017 at 12:30:40 -0500, Nicholas Jankowski wrote:
> > If you can help identify what Matlab does in this situation, and what
> > should be fixed and where, then maybe something can be done about it.
>
> Running the same script in Matlab 2016b produces no errors and the
> following output.
> Note Octave output from cheby2 is a row vector, Matlab is a column. I
> don't know if this matters.

Yes, that is https://savannah.gnu.org/bugs/?49318

> Matlab output follows. Octave 4.2.0,
> signal pkg v.1.3.2 ouput shown below that.
[…]
>   0.999999629730721 + 0.000860545419704i
>   0.999999629730721 - 0.000860545419704i
[…]
>   0.999999629730722 - 0.000860545419713i
>   0.999999629730722 + 0.000860545419685i

Thanks for running that. So yeah, looks like the bilinear function could
use some improvement here. If this last conjugate pair were correct,
there would be no error in cplxpair.

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank
On Thu, Jan 5, 2017 at 12:45 PM, Mike Miller <[hidden email]> wrote:

> On Thu, Jan 05, 2017 at 12:30:40 -0500, Nicholas Jankowski wrote:
>> > If you can help identify what Matlab does in this situation, and what
>> > should be fixed and where, then maybe something can be done about it.
>>
>> Running the same script in Matlab 2016b produces no errors and the
>> following output.
>> Note Octave output from cheby2 is a row vector, Matlab is a column. I
>> don't know if this matters.
>
> Yes, that is https://savannah.gnu.org/bugs/?49318
>
>> Matlab output follows. Octave 4.2.0,
>> signal pkg v.1.3.2 ouput shown below that.
> […]
>>   0.999999629730721 + 0.000860545419704i
>>   0.999999629730721 - 0.000860545419704i
> […]
>>   0.999999629730722 - 0.000860545419713i
>>   0.999999629730722 + 0.000860545419685i
>
> Thanks for running that. So yeah, looks like the bilinear function could
> use some improvement here. If this last conjugate pair were correct,
> there would be no error in cplxpair.
>
> --
> mike


Yes, Matlab's cplxpair errors the same way on Octave's z values as
well, so octave's cplxpair seems okay.

nickj

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

John W. Eaton-3
In reply to this post by hale812
On 01/05/2017 12:30 AM, hale812 wrote:

> This is really an awful kind of problem. Octave is becoming popular among
> engineers, but such fundamental problems, like failing signal toolbox, or
> crashing when computing 32bit numbers in long loops, reduce it's reputation
> to nothing. I just can not finish my work without rewriting everything in
> Fortran or C.
>
> 5 years guys! The problem is there longer than 5 years and nobody done a
> thing about it.

If you haven't already, I encourage you to read

   http://www.gnu.org/software/octave/support-expectations.html

jwe


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Mike Miller-4
In reply to this post by Mike Miller-4
On Thu, Jan 05, 2017 at 09:45:29 -0800, Mike Miller wrote:
> Thanks for running that. So yeah, looks like the bilinear function could
> use some improvement here. If this last conjugate pair were correct,
> there would be no error in cplxpair.

To clarify, to the OP in particular, this does not necessarily mean that
any work is now automatically going to be done to fix this.

If you care about this bug in particular, which it seems like you do,
you should expect to do some work to help get it fixed. If you don't at
the very least file a bug report, and preferably do some investigation
to help fix this, nothing will be fixed.

This thread will be forgotten by the end of this week or next. Someone
may very well come along in 2022 and ask why this is still not working
correctly if you don't do something now to help fix it.

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Ozzy Lash
I'm not sure that it is a problem in bilinear, it may be in sftrans.  If you do the pre-warping that is done in the cheby2 function for the digital case, and then run cheby2 with the "s" flag, the zeroes that result also don't pass the cplxpair test, unless you raise the tolerance to about 400*eps

>> pkg load signal
>> format long
>> stages=6;
>> reject=40;
>> fsamp = 1000;
>> fpL=0.52915085925565918945068233369966;
>> fpH=491.55617297547132693580351769924;
>> Ww=2/2*tan(pi*[fpL.*2./fsamp,fpH.*2./fsamp]/2)
Ww =

   1.66237798340326e-03   3.76885052761635e+01

>> [zw,pw,kw]=cheby2(stages,reject,Ww,"s");
>> [~,idx] = sort(abs(zw));
>> zw(idx).'
ans =

     0.000000000000000 -   0.000430272789501i
     0.000000000000000 +   0.000430272789515i
     0.000000000000000 +   0.001175504669177i
     0.000000000000000 -   0.001175504669177i
     0.000000000000000 +   0.001605738571502i
     0.000000000000000 -   0.001605738571502i
     0.000000000000000 +  39.017896505914251i
    -0.000000000000000 -  39.017896505914251i
     0.000000000000000 +  53.298419854299617i
    -0.000000000000000 -  53.298419854299624i
     0.000000000000000 + 145.611209739393132i
    -0.000000000000000 - 145.611209739393217i

>> cplxpair(zw)
error: cplxpair: could not pair all complex numbers
error: called from
    cplxpair at line 128 column 9


I'd guess that the issue may be in sftrans around this area:

        b = (C*(Fh-Fl)/2)./Sz;
        Sz = [b+sqrt(b.^2-Fh*Fl), b-sqrt(b.^2-Fh*Fl)];

when you have larrge and small magnitude values of zeroes, you probably lose precision.

Bill

On Thu, Jan 5, 2017 at 12:38 PM, Mike Miller <[hidden email]> wrote:
On Thu, Jan 05, 2017 at 09:45:29 -0800, Mike Miller wrote:
> Thanks for running that. So yeah, looks like the bilinear function could
> use some improvement here. If this last conjugate pair were correct,
> there would be no error in cplxpair.

To clarify, to the OP in particular, this does not necessarily mean that
any work is now automatically going to be done to fix this.

If you care about this bug in particular, which it seems like you do,
you should expect to do some work to help get it fixed. If you don't at
the very least file a bug report, and preferably do some investigation
to help fix this, nothing will be fixed.

This thread will be forgotten by the end of this week or next. Someone
may very well come along in 2022 and ask why this is still not working
correctly if you don't do something now to help fix it.

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

hale812
This post was updated on .
In reply to this post by John W. Eaton-3
Dear John.
If you haven't already, I encourage you to read
"https://www.gnu.org/software/octave/"
- Drop-in compatible with many Matlab scripts
- The Octave syntax is largely compatible with Matlab.
- packages for GNU Octave, similar to Matlab's toolboxes

So, this is a critical bug that should be fixed ASAP, while anyone intends to keep compatibility.
At least, the compatibility with common sense.

You know, the good code is the code that works unattended in specified constraints.

And in this case even working constraints are not specified. The error occurs in random and is not handled by any means. (Well, the random is weighted, but still it is random)
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Mike Miller-4
On Thu, Jan 05, 2017 at 15:55:00 -0800, hale812 wrote:

> Dear John.
> If you haven't already, I encourage you to read
> "https://www.gnu.org/software/octave/"
> - Drop-in compatible with many Matlab scripts
> - The Octave syntax is largely compatible with Matlab.
> - packages for GNU Octave, similar to Matlab's toolboxes
>
> So, this is a critical bug that should be fixed ASAP, while anyone intends
> to keep compatibility.
> At least, compatibility with common sense.

Ok. In case you're genuinely not aware, Octave is a volunteer project.
This is a discussion forum, not a support ticket. You are not a
customer. No one is obligated to fix anything. You are a participant in
a community, and if you care about this particular bug, you will make an
effort to try and help fix it. We hope that you will.

Cheers,

--
mike

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank
In reply to this post by hale812
On Thu, Jan 5, 2017 at 6:55 PM, hale812 <[hidden email]> wrote:
> Dear John.
> If you haven't already, I encourage you to read
> "https://www.gnu.org/software/octave/"
> - Drop-in compatible with many Matlab scripts
> - The Octave syntax is largely compatible with Matlab.
> - packages for GNU Octave, similar to Matlab's toolboxes
>
> So, this is a critical bug that should be fixed ASAP

We do have bug reports going back to 2011, I've recently chimed in on
a few to try to bring some of the older ones to a close. I must be
missing the bug report you helped log the first time with your
detailed knowledge of the error so this wouldn't fall off the radar.
Can you point me to it? it should be over at http://bugs.octave.org

If you can't find it (the search is less than optimal) maybe you'd
consider starting a new one with a minimal working example [1] [2] we
can use to gradually work through debugging the error? Shouldn't take
more time than you've already spent on the help emails.

[1] https://en.wikipedia.org/wiki/Minimal_Working_Example
[2] http://stackoverflow.com/help/mcve

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

Ozzy Lash
In reply to this post by Ozzy Lash


On Thu, Jan 5, 2017 at 1:30 PM, Ozzy Lash <[hidden email]> wrote:
I'm not sure that it is a problem in bilinear, it may be in sftrans.  If you do the pre-warping that is done in the cheby2 function for the digital case, and then run cheby2 with the "s" flag, the zeroes that result also don't pass the cplxpair test, unless you raise the tolerance to about 400*eps


Looking a bit deeper, I think that the issue may be in the way the zeros of the nth order chebychev type 2 filter are being calculated.  The zeros of the prototype, while they successfully can be paired with cplxpair, they are not complex conjugates.

For the case of an even order, the zeros are calculated pretty simply as:

theta = pi*([1:n]-0.5)/n

z=1i./cos(theta)

In the odd order case, the zero at pi/2 is omitted, so it is a little more complex.

for n=6, you get:

>> n=6
n =  6
>> theta=pi*([1:n]-0.5)/n
theta =

 Columns 1 through 4:

   0.261799387799149   0.785398163397448   1.308996938995747   1.832595714594046

 Columns 5 and 6:

   2.356194490192345   2.879793265790644

>> z=1i./cos(theta)
z =

 Columns 1 and 2:

   0.00000000000000 + 1.03527618041008i   0.00000000000000 + 1.41421356237309i

 Columns 3 and 4:

   0.00000000000000 + 3.86370330515627i  -0.00000000000000 - 3.86370330515628i

 Columns 5 and 6:

  -0.00000000000000 - 1.41421356237310i  -0.00000000000000 - 1.03527618041008i

>> cplxpair(z)
ans =

 Columns 1 and 2:

  -0.00000000000000 - 1.03527618041008i   0.00000000000000 + 1.03527618041008i

 Columns 3 and 4:

  -0.00000000000000 - 3.86370330515628i   0.00000000000000 + 3.86370330515627i

 Columns 5 and 6:

  -0.00000000000000 - 1.41421356237310i   0.00000000000000 + 1.41421356237309i


instead of calculating all of the zeros, I wonder if it would make sens to calculate the zeros for values of theta from 0 up to but not including pi/2, and then conjugating those zeros, so that all of the zeros are complex conjugate pairs. 

Here is a patch that makes this change. 

--- /home/lash/octave/signal-1.3.2/cheby2.m    2017-01-05 22:10:48.630452347 -0600
+++ cheby2.m    2017-01-05 22:37:24.415014904 -0600
@@ -113,10 +113,11 @@
   beta = cosh (phi) * cos (theta);
   if (rem (n, 2))
     ## drop theta==pi/2 since it results in a zero at infinity
-    zero = 1i * C ./ cos (theta([1:(n - 1) / 2, (n + 3) / 2:n]));
+    zero = 1i * C ./ cos (theta([1:(n - 1) / 2]));
   else
-    zero = 1i * C ./ cos (theta);
+    zero = 1i * C ./ cos (theta([1:n/2]));
   endif
+  zero = [zero, conj(zero)];
   pole = C ./ (alpha.^2 + beta.^2) .* (alpha - 1i * beta);
 
   ## Compensate for amplitude at s=0


Is there a bug report for this issue?  I can add this to the bug report if so.

I'm not sure that this is the best solution, and there may be issues with the poles as well as the zeros.

Bill

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

hale812
No, there is no bug report. I was using Scilab and Matlab before, but at the moment moved to Octave. So, I do not have an access to bug-tracking channel yet.
Thank you very much for looking deep into the problem!
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

hale812
In reply to this post by Mike Miller-4
Mike. I can wait patiently...
But 5 years Mike! More than 5 years!

I would like to start the bug ticket again, but as I told it, I have no access to bug tracking yet. (I've sent a request)

And I told you earlier, it is not a matter of excuses. (since a lazy excuse, that is)

It is a matter of common sense. If you release a good code (not alpha and not beta), it should work in common constraints. Or the constraints should be defined, so less educated people (like me) could track the values and make workarounds.

But you see, some computation scripts can grow to megabytes of perfectly correct code, worse, if there is I/O, GUI and multithreading.
And with such fundamental bugs... Octave is just gets devalued.
Who needs all your great and voluntary work, if there's a bomb buried under it?
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank

On Jan 6, 2017 3:20 AM, "hale812" <[hidden email]> wrote:
>
> I would like to start the bug ticket again, but as I told it, I have no
> access to bug tracking yet. (I've sent a request)

Everyone has access to the bug tracker.

http://bugs.octave.org

Anyone can submit a bug report. A login account is optional. Give it a try. Along the top right of the page there's a Bugs menu with Submit New as an option. Fill in what you're able and the process starts.

Also, this is a hard stop error that has so far been brought up 2x in 5 years. It also is in a package, not core octave. This may seem critical to you, but if it was a roadblock to more people odds are that it would have been heard about more often.

Nick j.


_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank
In reply to this post by Ozzy Lash
On Thu, Jan 5, 2017 at 11:44 PM, Ozzy Lash <[hidden email]> wrote:

> Here is a patch that makes this change.
>
> --- /home/lash/octave/signal-1.3.2/cheby2.m    2017-01-05 22:10:48.630452347
> -0600
> +++ cheby2.m    2017-01-05 22:37:24.415014904 -0600
> @@ -113,10 +113,11 @@
>    beta = cosh (phi) * cos (theta);
>    if (rem (n, 2))
>      ## drop theta==pi/2 since it results in a zero at infinity
> -    zero = 1i * C ./ cos (theta([1:(n - 1) / 2, (n + 3) / 2:n]));
> +    zero = 1i * C ./ cos (theta([1:(n - 1) / 2]));
>    else
> -    zero = 1i * C ./ cos (theta);
> +    zero = 1i * C ./ cos (theta([1:n/2]));
>    endif
> +  zero = [zero, conj(zero)];
>    pole = C ./ (alpha.^2 + beta.^2) .* (alpha - 1i * beta);
>
>    ## Compensate for amplitude at s=0
>
>

This does seem to fix the issue, but I also can't tell whether it's a
'correct' fix:


-----------------------------------
>> [z,p,k]=cheby2(6,40,[fpL.*2 ./1000,fpH.*2./1000]);z'
ans =

  -0.998687143937134 - 0.051224881988043i
  -0.999296201278521 - 0.037511359750325i
  -0.999905676503125 - 0.013734558486826i
   0.999994843220576 - 0.003211468862565i
   0.999997236381364 - 0.002351006089708i
   0.999999629730722 - 0.000860545419685i
   0.999994843220576 + 0.003211468862565i
   0.999997236381364 + 0.002351006089708i
   0.999999629730722 + 0.000860545419685i
  -0.998687143937134 + 0.051224881988043i
  -0.999296201278521 + 0.037511359750325i
  -0.999905676503125 + 0.013734558486826i

>> cplxpair(z)'
ans =

  -0.999905676503125 + 0.013734558486826i
  -0.999905676503125 - 0.013734558486826i
  -0.999296201278521 + 0.037511359750325i
  -0.999296201278521 - 0.037511359750325i
  -0.998687143937134 + 0.051224881988043i
  -0.998687143937134 - 0.051224881988043i
   0.999994843220576 + 0.003211468862565i
   0.999994843220576 - 0.003211468862565i
   0.999997236381364 + 0.002351006089708i
   0.999997236381364 - 0.002351006089708i
   0.999999629730722 + 0.000860545419685i
   0.999999629730722 - 0.000860545419685i

>> [sos, g] = zp2sos (z, p, k)
sos =

 Columns 1 through 3:

   1.000000000000000   1.999811353006250   1.000000000000000
   1.000000000000000   1.998592402557041   1.000000000000000
   1.000000000000000   1.997374287874269   1.000000000000000
   1.000000000000000  -1.999989686441153   1.000000000000000
   1.000000000000000  -1.999994472762729   1.000000000000000
   1.000000000000000  -1.999999259461444   1.000000000000000

 Columns 4 through 6:

   1.000000000000000   1.967506243219743   0.972881387742366
   1.000000000000000   1.923500197249973   0.927580075649338
   1.000000000000000   1.899456851014250   0.902325817320743
   1.000000000000000  -1.993570273311973   0.993582084086452
   1.000000000000000  -1.995281755784743   0.995298350815642
   1.000000000000000  -1.998255341271945   0.998276741178970

g =  0.896584733013340

-------------------------

recalling from earlier:

Matlab 2016b:
>>  [sos, g] = zp2sos (z, p, k)
sos =
  Columns 1 through 3
   1.000000000000000   1.999811353006249   0.999999999999999
   1.000000000000000   1.998592402557042   1.000000000000000
   1.000000000000000   1.997374287874268   0.999999999999998
   1.000000000000000  -1.999999259461443   0.999999999999999
   1.000000000000000  -1.999994472762730   1.000000000000001
   1.000000000000000  -1.999989686441154   1.000000000000001

  Columns 4 through 6
  1.000000000000000   1.899456851014267   0.902325817320760
  1.000000000000000   1.923500197249950   0.927580075649315
  1.000000000000000   1.967506243219755   0.972881387742377
  1.000000000000000  -1.993570273311965   0.993582084086443
  1.000000000000000  -1.995281755784758   0.995298350815657
  1.000000000000000  -1.998255341271938   0.998276741178964
g =
   0.896584733013338

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: zp2sos error: cplxpair: could not pair all complex numbers

nrjank
again, no intuitive sense of what I'm looking at, but how important is
the sort order for sos? The values all seem the same out to 1E-14 or
so, but matlab seems to have them sorted differently.

_______________________________________________
Help-octave mailing list
[hidden email]
https://lists.gnu.org/mailman/listinfo/help-octave
12