firls.m, part 2

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
32 messages Options
12
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

firls.m, part 2

je suis
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

NJank
On Sun, Jun 18, 2017 at 10:55 AM, je suis <[hidden email]> wrote:
> 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

Well, a number of people are working on a number of their own
projects, it appears that you were point person on this one.  Looking
back through your last request for input was:

"Please check the two scripts for whatever input you can think of."

Well, this isn't a routine that I currently have an immediate
technical use for, so the only inputs I'd have to test are the
examples we've already gone through. With the modified expint, do they
reproduce well? That was an interesting find and might warrant its own
discussion. 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.

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.

I'm thinking the best thing might now be to create two separate
submissions to the Octave bug tracker, one for firls (if you haven't
done this already) and one for expint at bugs.octave.org. You can link
the github archive and mention the dependency between the two bug
reports.

 I've cc'd Mike, who is maintaining the signals package.  He may have
suggestions on what else is needed for this function prior to
inclusion. There is no single maintainer for the specfun package,
which contains expint, but maybe he has suggestions there.

Nick J.

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

Re: firls.m, part 2

Mike Miller-4
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
In reply to this post by NJank
> Well, a number of people are working on a number of their own
> projects, it appears that you were point person on this one.

I'm not sure how to translate this: "you were point person on this one".

> Well, this isn't a routine that I currently have an immediate
> technical use for, so the only inputs I'd have to test are the
> examples we've already gone through. With the modified expint, do they
> reproduce well? That was an interesting find and might warrant its own
> discussion. 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.
>
> 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.

You're right, I'm sorry, I completely forgot about those. I'll modify
the script(s).

> I'm thinking the best thing might now be to create two separate
> submissions to the Octave bug tracker, one for firls (if you haven't
> done this already) and one for expint at bugs.octave.org. You can link
> the github archive and mention the dependency between the two bug
> reports.

Then I disregarded the etiquette, for which I apologize.

>  I've cc'd Mike, who is maintaining the signals package.  He may have
> suggestions on what else is needed for this function prior to
> inclusion. There is no single maintainer for the specfun package,
> which contains expint, but maybe he has suggestions there.

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.

Hope to hear good news soon,
Vlad

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

Re: firls.m, part 2

Mike Miller-4
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
> 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

Dmitri A. Sergatskov
On Mon, Jun 19, 2017 at 11:23 AM, je suis <[hidden email]> 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 <cmath>
#include <iostream>
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.c

Dmitri.
--





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

Re: firls.m, part 2

Mike Miller-4
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
> 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

Mike Miller-4
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
In reply to this post by Dmitri A. Sergatskov
> ​expint() is in C++17, so gcc7+ supports it:
>
> ​cat c_expint.cc
> #include <cmath>
> #include <iostream>
> 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.c

I 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

Francesco Potortì
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

NJank
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

NJank
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
@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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

Mike Miller-4
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

NJank
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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: firls.m, part 2

je suis
In reply to this post by NJank
Using the link for expint() shows spot on results, compared to my
cooked up version.

> 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.

There are some checks in the beginning, those should take care of too
few, or too many arguments, correct string arguments, and correct
numeric ones (thought not for N). Should I delete those and replace
them with the %! tyoe checks you gave as example?

Vlad

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