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:
and the cplxpair_eqlz.m is
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.
|
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
But the values that SciPy computes for the same filter design call behave when pasted into Octave:
|
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);
|
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:
dunno if there's a 'nice round number' set of bug-tripping values we can use as with the butter example.
|
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.
|
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.
|
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.
|
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.
|
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) ...
|
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 _?_
|
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:
in Matlab 2016b, the output produced is:
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?
|
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)
|
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.
|
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.
|
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.
|