ENH: signal: add `invfreqs` and `invfreqz`

[IADSJohn]

Is anyone still working on an invfreqs() funtion for Scipy? I saw a close issue ticket on the matter...We could really use one.

[hgomersall]

Also, invfreqz!

[IADSJohn]

New function?

John Bretz

IADS Development & Support

(661)277-5410

(661)273-7003 x205

From: Henry Gomersall [mailto:notifications@github.com] Sent: Tuesday, February 02, 2016 8:00 AM To: scipy/scipy scipy@noreply.github.com Cc: IADSJohn jcbretz@symvionics.com Subject: Re: [scipy] invfreqs() Function Needed... (#4497)

Also, invfreqz!

— Reply to this email directly or view it on GitHub https://github.com/scipy/scipy/issues/4497#issuecomment-178653526 . https://github.com/notifications/beacon/AKXgfcJLrwBCBoC8g1cKAOKYIrKKgmxXks5pgMoLgaJpZM4DcYUI.gif

[rgommers]

For reference, that was gh-920 (closed issue because patch had license issues). No one is working on this as far as I know, but both invfreqs and invfreqz would be good to have. Contributions very welcome.

[hgomersall]

At the minute I'm struggling to find any decent references for this. Any anybody knows of would be much appreciated and would likely lead to a patch.

The Matlab implementation has stability issues (numerical whilst running, as well as the derived filters) or possibly documentation problems. It's also based on very old theory.

[endolith]

Did you compare with the Octave version, which does both? http://octave.sourceforge.net/signal/function/invfreq.html

[hgomersall]

Thanks for that, I'd assumed it was the same algorithm initially, but it seems actually to perform differently. I'll have a play with it later on today. I don't think is has any stability enforcing constraints, but I might be wrong. Perhaps I can combine what I learn there with some more recent papers...

[abhijeet-pandey]

Any activity on adding - invfreqs and invfreqz to scipy.signal?

[tylerjereddy]

@abhijeet-pandey doesn't look like it--I did a search of open PRs for both terms and nothing came up. Plus there's no cross-linked PR above, and git grep comes up negative for both terms, so it seems like the situation is still that this would be a welcome PR/addition.

[josmithiii]

Hi, I'm the original author of invfreq*, and I hereby happily upgrade license to my contributions to whatever scipy is using, or MIT license, etc. I wrote the original version as a favor for my friend John Little, founder of The Mathworks, and did not get paid for it, so I still own the copyright to what I wrote. I also started the octave-forge version. I expect anyone who contributed to that version will go along with this license upgrade as well. The algorithm is described here.

[lucascolley]

Thanks @josmithiii !

Would anyone like to send a PR?

[josmithiii]

I'm happy to help with testing and debugging. I found this discussion because I need it in Python myself!

[abhijeet-pandey]

I need a working example code that implements the invfreq algorithm and I can try to port it to Scipy. @josmithiii - Thanks for sharing the resource. I am unwilling to directly use algorithm at https://ccrma.stanford.edu/~jos/filters/FFT_Based_Equation_Error_Method.html to code in Scipy as I see a number of special cases, in Octave, depending on count of variables, and am not confident about coding it all out, using just the algorithm page as reference, without introducing any significant bug. I am already having trouble interpreting use of array slicing with negative indices. A slice -1:1 in a length 5 vector - does it produce [-1, 0, 1] or [-1, -2, -3, -4, -5] since index -5 is the same as index 1? There is a slight typo Ruy = Ryu.T.

In short, I would need more lower level details. I had a look at the Octave implementation of the invfreq algorithm to see if it can be directly ported over to Scipy. However, I am not convinced that it has been implemented correctly. For e.g. I was unable to reproduce the results at https://au.mathworks.com/help/signal/ref/invfreqz.html in Octave. Relevant to this issue, I see there was an open bug reported at https://savannah.gnu.org/bugs/?52604. I understand that because of differences in algorithms (between Matlab and Octave), the results from two codes can be different but simultaneously correct. However, the order of magnitude of b coefficients returned by invfreq by Octave (for the Matlab example above, which was of the order of 10^13) made me suspicious of the whole implementation.

Secondly, I looked at another implementation of the algorithm, this time in python - https://github.com/yurochka/dsptools/blob/master/dsptools/invfreqz.py and I can't reproduce the original Matlab results either. If you could provide any working example implementation of invfreq in any language of your choice, I will be happy to port that to Scipy.

[josmithiii]

Thanks for the rapid look!

I would prefer that we start with the Octave code. What I sent was tested to a decent extent, but I see that you've found a failing case, so I need to look at that. I assume you are just using Octave with -forge invfreqz installed. I am unable to do that myself because the control package fails to link on my Mac M1 (ld: library 'octinterp' not found), and it's apparently not because I didn't install it, but an incorrect build recipe. I'll try to mdfind the source and copy it to a new name.

Huge results are normally due to a singular matrix inversion. There are also differences of course when the algorithm is not normalized to be the same.

The open GNU issue is a "missing feature" complaint. The iteration feature is what is called the Steiglitz McBride algorithm, as my write-up mentions. I've never implemented it. Are scipy functions allowed to implement only some of the expected features or is that prohibited? I know Octave allows it. It's not hard to add SM, but it introduces a lot of UI complexity because it's iterative, and what if it doesn't converge, what if the starting poles are unstable, etc. The base invfreq* algorithm is just a bone-simple least-squares solution - one iteration of Newton's method. The only way it can get into trouble is to have a singular matrix inverse on its hands (which does happen pretty often, especially for challenging filter design problems). It also does not guarantee filter stability. I have a blog on this issue and how to hopefully avoid it. As always, the more you know, the more effective are your tools. I'm happy to incorporate that info to the scipy doc if someone points me to how.

If a vector V is length 5, and indexing is zero-based (unlike matlab), then V[-1:1] means ( V[4], V[0], V[1] ), so I guess that's what you're describing as [-1, 0, 1].

I'll also check out the dsptools version, but I lean toward fixing and enhancing versions based on my starting code.

Thanks for reporting the typo - fixed in the source. I also fixed that final ','.

I'll let you know when I know more. Also let me know if I can access the tests you've done already, as it would probably save me time (always the most precious resource, especially at my age).

[josmithiii]

I got Octave's invfreqz going by abandoning brew and installing octave-signal using MacPorts - smooth sailing. The Mathworks example is garbage, in my opinion, because the filter poles are unstable, and worse, two are right on the unit circle:

octave:1> format long
octave:2> a = [1   2   3   2   1   4];
octave:3> abs(roots(a))
ans =

   1.453397651516403
   1.658967081916993
   1.658967081916993
   1.000000000000001
   1.000000000000001

Therefore, the output of freqz does not even exist, to numerical precision. I have no interest in trying to match a degenerate and unmotivated differing example. There are problems enough trying to do things that are worth doing. I like how Octave warns about ill conditioning. I think that's all we need on the numerical failure front.

So, I now have /opt/local/share/octave/packages/signal-1.4.5/invfreqz.m and I can improve it if we find any real need to do so (or a new failing example that's worth something).

From the source comments, it looks like only Andrew Fitting has added code beyond what I wrote. Should we try to reach him regarding the license upgrade before proceeding? How about we have Claude 3.5 Sonnet translate it to Python for us? I know this is controversial, but according to my limited understanding of copyright law, that should be technically allowed even if we can't find Andrew (who I do expect not to object anyway). Fun times . . .

[abhijeet-pandey]

inverseFreq.zip Sorry, I don't have much knowledge on the license issue. Attached are the Octave/Matlab example problem and the Octave algorithm implemented in python which produces the same coefficients except for a few coefficients that are may be affected by inversion of badly conditioned matrix. It is incomplete but should be a good starting point. Hope it is of some use. I would really like to get this to a point that it is able to replicate results of the Matlab example. By the way, I am a huge admirer of https://ccrma.stanford.edu/~jos/filters/ - which I think is the best exposition on digital filters! Thank you.

[lucascolley]

From the source comments, it looks like only Andrew Fitting has added code beyond what I wrote. Should we try to reach him regarding the license upgrade before proceeding?

That sounds sensible.

[ilayn]

@josmithiii @abhijeet-pandey Thank you for addressing this issue once again.

Please make sure that you are not drawing, even being inspired by, other GPL licensed code. It is unfortunately a real issue. Hence we have to be very careful with Octave code, we can't keep investingating that code. Here though you have originated the implementation it lives in that code base and further modified already. Thus it does not still belong to you.

Instead of working on licensed code, if you have the code again in matlab-lang or any other language, even pseudo-code, we can convert it to Python, Cython or C or whatever, that's not a problem I can help with that. We just need a clean, not necessarily performant, code to start with and we can squeeze the nuts and bolts afterwards. We just need to cut ties with GPL code. I am familiar with the subject matter hence we can work on an independent implementation.

[ilayn]

I checked the matlab help page for invfreqs and they are doing a direct OE based computation and hence should be doable (famous last words).

[josmithiii]

Thanks for the helpful and kind comments! Ok, I will abandon the Octave code and will try to help @abhijeet-pandey finish his fresh implementation. (I doubt I will find my original 1986 implementation, but I'll do some find searches.) I'm also curious to see what Claude 3.5 (and soon Orion) can do with my LaTeX description. Claude is very fluent in LaTeX, and can do basic math in terms of it.

So, to be clear, we can only start from the public API, right? How about the help info? Or doc info? Surely this has all been tested repeatedly in court by now. As I recall, Java was a big test case, but I didn't follow the full story. Also, AT&T failed to bring down GNU and Linux. It appears APIs are safe to reimplement, but I'm not sure about man pages, help, doc, etc.

Is there a written scipy policy document somewhere that spells all this out?

[josmithiii]

Regarding the difficult example, try changing the 64 to 63 or 66. LOL! There are multiple ways to address the issues [(1) poles on the unit circle, and (2) non-minimum-phase input spectrum], but the general category is "Garbage In, Garbage Out". We could put such mitigations in the code (and win on all three cases), but that goes outside the scope of invfreqz, and should probably have two new boolean arguments such as convert_to_minimum_phase_if_need_be and factor_out_large_spectral_peaks_and_treat_them_separately or equivalent.

[lucascolley]

Is there a written scipy policy document somewhere that spells all this out?

A relevant section from our docs is:

PRs are often submitted with content copied or derived from unlicensed code or code from a default license that is not compatible with SciPy’s license. For instance, code published on StackOverflow is covered by a CC-BY-SA license, which is not compatible due to the share-alike clause. These contributions cannot be accepted for inclusion in SciPy unless the original code author is willing to (re)license their code under the modified BSD (or compatible) license. If the original author agrees, add a comment saying so to the source files and forward the relevant communication to the scipy-dev forum.

Another common occurrence is for code to be translated or derived from code in R, Octave (both GPL-licensed) or a commercial application. Such code also cannot be included in SciPy. Simply implementing functionality with the same API as found in R/Octave/… is fine though, as long as the author doesn’t look at the original incompatibly-licensed source code.

That doesn't spell out everything in detail, but I think that copying parts of documentation would be similarly problematic.

[josmithiii]

Ok, here is my result from having Claude translate my online description of the algorithm (it did very well): invfreqz_jos.py.tgz It passes my first several tests, but needs more testing, and expansion to more argument cases. I think I should work on Steiglitz-McBride iterations next.

[josmithiii]

I created a scipy fork here. Branch joscontains some temporary files scipy/scipy/signal/*_jos.py from which a proposed invfreqz for _filter_design.py can be derived (after some obvious post-development cleanup). My initial tests pass for the non-iterative case, but more are needed. Please clone, create your own branch, and add whatever features/bug-fixes you feel would be needed or nice (especially @abhijeet-pandey who already started on some code).

I have ideas for improving the iterative case that I plan to try out soon (analogous to a learning-rate schedule in deep learning). I also have alternative methods to propose, such as two variations of Prony's method, so a new method parameter would be added. Its default would be 0 giving the usual frequency-domain equation-error method. Please let me know if I threaten to duplicate something already somewhere else in the Python universe, or if we don't want so many classic filter-design methods. I'm assuming _filter_design.py contains all that exists so far for IIR filter design.