# firls.m, part 2

32 messages
12
Open this post in threaded view
|
Report Content as Inappropriate

## firls.m, part 2

 Hello The previous thread got too long, so I thought I'd start a new one (sorry and let me know if it's the wrong way). I managed to make ...something, on github: https://github.com/archbugaboo/Octave_signal. The new script is modified and should perform better now, but there are still minor differences. I suppose they have their way to circumvent cos(0)/0 and Ci(0), I have mine? Please test it and let me know if you think it suits you. Vlad PS: I'll be using the code for a program of mine, it's only on my laptop, for now, but if it will be online, it will be open source (Qt). I hope there isn't any clash in licenses. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 I have discovered the source of the errors: expint(). For purely imaginary arguments, as their value increases, the results come more and more worse. I have made a separate version, E1(x), which uses code adapted from Numerical recipes and, not least, when run with this line: x=(1:0.1:100)*1i; tic; a=E1(x); toc, tic; b=expint(x); toc the timings are almost the same: 0.76s (E1) and 0.72 (expint). This for eps=1e-15. I gave uploaded E1.m in the same location on github and, after some more testing with the 3 examples I have with 1/f^2 weighting [mostly plotting diff(abs(fft(h)))) ],  I dare conclude that my numbers are better than Matlab's. Please check the two scripts for whatever input you can think of. When using 1/f^2 weighting, don't forget that there will be similar numbers, with larger differences towards the ends of the impulse, and with a dew common decimals towards the middle, but I dare you to plot the derivative of the FFT for both: you can clearly see the exponential rise in the ripple and the (almost) zero ripple close to DC, in my results. I'm sorry, I can only stay modest for so long. :-) Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 Given the answers in the previous threads I assumed there is some interest in this, but now that days seem to pass, I have to ask: is there anymore left? Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 On Mon, Jun 19, 2017 at 11:43:36 -0400, Nicholas Jankowski wrote: > Well, a number of people are working on a number of their own > projects, it appears that you were point person on this one. Yeah, very often the number of people interested in a function enough to work on fixing it or improving it is between zero and one. > If there actually is a deficiency in expint, and your > version improves it I think fixing the main one would be preferable to > just having your fix as a subfunction to firls. Absolutely. If there is a problem with Octave's expint for certain types of inputs, please consider fixing that there instead of introducing a workaround function. > I notice there are still no built-in tests at the end of the current > version.  I'd recommend adding a bunch of those for basic input/output > form checking as well as a few (simple if possible) expected numerical > outputs. I think one of my past emails may have included some > examples. Built-in tests are a must, firls.m has none and any changes must include some tests to validate that it's working as expected. > There is no single maintainer for the specfun package, > which contains expint, but maybe he has suggestions there. The expint function used to be in the specfun package, it is now part of Octave. -- mike _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 On Mon, Jun 19, 2017 at 16:07:55 +0000, je suis wrote: > Thank you very much for your help so far, I'll try to continue on the > right path. I don't know if you got to that part, but the expint() > implementation is adapted from Numerical Recipes, I hope that doesn't > come against any licenses. I suspect that it might, depending on what you mean by "adapted from". You may want to take down your code wherever it is hosted publicly. Take a look at http://numerical.recipes/licenses/redistribute.html. -- mike _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 > I suspect that it might, depending on what you mean by "adapted from". > You may want to take down your code wherever it is hosted publicly. Take > a look at http://numerical.recipes/licenses/redistribute.html. I see, but the underlying mathematics is based on Abramowitx & Stegun. If I modify the code, it should be safe then, shouldn't it? Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 On Mon, Jun 19, 2017 at 11:23 AM, je suis wrote:> I suspect that it might, depending on what you mean by "adapted from". > You may want to take down your code wherever it is hosted publicly. Take > a look at http://numerical.recipes/licenses/redistribute.html. I see, but the underlying mathematics is based on Abramowitx & Stegun. If I modify the code, it should be safe then, shouldn't it? Vlad ​expint() is in C++17, so gcc7+ supports it:​cat c_expint.cc #include #include int main(){    std::cout << "Ei(0) = " << std::expint(0) << '\n'              << "Ei(1) = " << std::expint(1) << '\n'              << "Gompetz constant = " << -std::exp(1)*std::expint(-1) << '\n';}g++ -std=c++1z c_expint.cc./a.out Ei(0) = -infEi(1) = 1.89512Gompetz constant = 0.596347Also there is a free implementation in e.g. GSL:https://github.com/ampl/gsl/blob/master/specfunc/expint.cDmitri.-- _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 In reply to this post by je suis On Mon, Jun 19, 2017 at 16:23:45 +0000, je suis wrote: > > I suspect that it might, depending on what you mean by "adapted from". > > You may want to take down your code wherever it is hosted publicly. Take > > a look at http://numerical.recipes/licenses/redistribute.html. > > I see, but the underlying mathematics is based on Abramowitx & Stegun. > If I modify the code, it should be safe then, shouldn't it? There is no right answer and I am not a lawyer. I also haven't read NR, I haven't looked at the relevant A&S section, and I haven't looked at your code. The expint function in Octave has comments indicating parts of it are also based on A&S. My recommendation would be to start from the current expint function and make the necessary changes to it to fix the problem you have found, based on what you have learned about the algorithm. -- mike _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 > There is no right answer and I am not a lawyer. I also haven't read NR, > I haven't looked at the relevant A&S section, and I haven't looked at > your code. The expint function in Octave has comments indicating parts > of it are also based on A&S. My recommendation would be to start from > the current expint function and make the necessary changes to it to fix > the problem you have found, based on what you have learned about the > algorithm. Well, if you say your own code has parts based on A&S, then it should be safe (as I imagine it is for so many others that use A&S as a reference, or mathematics, for that matter). Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 In reply to this post by Mike Miller-4 On Mon, Jun 19, 2017 at 09:35:35 -0700, Mike Miller wrote: > the current expint function and make the necessary changes to it to fix > the problem you have found, based on what you have learned about the > algorithm. Also make sure you have actually tested the correct current version of the expint function. I haven't seen what version of Octave you are using, but there have been three different rewrites of expint over time. The most recent revision was completely rewritten in January 2017 and is not part of any Octave release yet. Please test your cases against it to see whether the problem you encountered has already been fixed. -- mike _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 In reply to this post by Dmitri A. Sergatskov > ​expint() is in C++17, so gcc7+ supports it: > > ​cat c_expint.cc > #include > #include > int main() > { >     std::cout << "Ei(0) = " << std::expint(0) << '\n' >               << "Ei(1) = " << std::expint(1) << '\n' >               << "Gompetz constant = " << -std::exp(1)*std::expint(-1) << > '\n'; > } > > g++ -std=c++1z c_expint.cc > ./a.out > Ei(0) = -inf > Ei(1) = 1.89512 > Gompetz constant = 0.596347 > > Also there is a free implementation in e.g. GSL: > https://github.com/ampl/gsl/blob/master/specfunc/expint.cI didn't dare use (or see about) C++17 (but I am relatively new with this). Thank you for letting me know of this. Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 In reply to this post by je suis Mike Miller: >> I suspect that it might, depending on what you mean by "adapted from". >> You may want to take down your code wherever it is hosted publicly. Take >> a look at http://numerical.recipes/licenses/redistribute.html. Vlad: >I see, but the underlying mathematics is based on Abramowitx & Stegun. >If I modify the code, it should be safe then, shouldn't it? The problem is not with the mathematics, is with the code.  Copyright law applies to the code: if you copy it, translate it or modify it, you are using copyrigted code, and you should have a licence to do so.  The above page is concerned with the case when you use their code. If you follow a math concept, a description, a conceptual flow, and then implement it yourself, you have no problems: you own the copyright to the code you write andyou can license it as you wish. -- Francesco Potortì (ricercatore)        Voice:  +39.050.621.3058 ISTI - Area della ricerca CNR          Mobile: +39.348.8283.107 via G. Moruzzi 1, I-56124 Pisa         Skype:  wnlabisti (entrance 20, 1st floor, room C71)     Web:    http://fly.isti.cnr.it_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 In reply to this post by Mike Miller-4 On Mon, Jun 19, 2017 at 12:44 PM, Mike Miller <[hidden email]> wrote: > On Mon, Jun 19, 2017 at 09:35:35 -0700, Mike Miller wrote: > The most recent revision was completely rewritten in January 2017 and is > not part of any Octave release yet. Please test your cases against it to > see whether the problem you encountered has already been fixed. Vlad, here's a link [1] to the more recent expint. you can just save it in your working folder or overwrite the one in Octave. the original bug report [2] also called out large errors with imaginary inputs. If this also fixes the problem I would suggest just making use of it. That'll avoid any copyright concerns with Numerical Recipes. [1] http://hg.savannah.gnu.org/hgweb/octave/file/c56d30ea6cd4/scripts/specfun/expint.m[2] https://savannah.gnu.org/bugs/?47738_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 On Mon, Jun 19, 2017 at 3:12 PM, Nicholas Jankowski <[hidden email]> wrote: > On Mon, Jun 19, 2017 at 12:44 PM, Mike Miller <[hidden email]> wrote: >> On Mon, Jun 19, 2017 at 09:35:35 -0700, Mike Miller wrote: > > Vlad, here's a link [1] to the more recent expint. you can just save a better link, this will give you the most recent version if any changes have been made since the update back in Jan. http://hg.savannah.gnu.org/hgweb/octave/file/tip/scripts/specfun/expint.m_______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 @all The day starts busy for me, but by tonight I should have done most, if not all. I've skimmed over the linked expint.m, but I see the same Numerical Recipes influence (even if it is Fortran). If all it takes is a modification of the code to not make it look like it has been adapted, then it should be fine, since, how many ways can you reinterpret A&S formulas? If not, the first one that makes a program out of A&S holds dominion over everyone. Maybe I'm too convoluted right now. Also, if I find a way to improve the linked expint.m script (unlikely, but you never know), is it alright to do so? Vlad _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 On Tue, Jun 20, 2017 at 05:03:39 +0000, je suis wrote: > The day starts busy for me, but by tonight I should have done most, if > not all. I've skimmed over the linked expint.m, but I see the same > Numerical Recipes influence (even if it is Fortran). If all it takes > is a modification of the code to not make it look like it has been > adapted, then it should be fine, since, how many ways can you > reinterpret A&S formulas? If not, the first one that makes a program > out of A&S holds dominion over everyone. Maybe I'm too convoluted > right now. Also, if I find a way to improve the linked expint.m script > (unlikely, but you never know), is it alright to do so? My main point was that expint has already been updated in Octave, and you have been working with an older revision of the function. So please test your firls function against the current expint to see if it works for you. If there are any remaining problems with the current expint, yes of course, we would be happy if you report a bug and help fix it if you can. -- mike _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave
Open this post in threaded view
|
Report Content as Inappropriate

## Re: firls.m, part 2

 On Tue, Jun 20, 2017 at 11:39 AM, Mike Miller <[hidden email]> wrote: > On Tue, Jun 20, 2017 at 05:03:39 +0000, je suis wrote: > So please test your firls function against the current expint to see if it works for you. I would add that you should take a peek at the test functions at the end of the newer version of expint. It has good examples testing both the basic input/output forms as well as some expected numerical output to a certain level of precision. _______________________________________________ Help-octave mailing list [hidden email] https://lists.gnu.org/mailman/listinfo/help-octave