John Ridges
2009-Jul-22 20:49 UTC
[Speex-dev] A technical question about the speex preprocessor.
By my reckoning the confluent hypergoemetric functions should have the following values: M(-.25;1;-.5) = 1.11433 M(-.25;1;-1) = 1.21088 M(-.25;1;-1.5) = 1.29385 M(-.25;1;-2) = 1.36627 M(-.25;1;-2.5) = 1.43038 M(-.25;1;-3) = 1.48784 M(-.25;1;-3.5) = 1.53988 M(-.25;1;-4) = 1.58747 M(-.25;1;-4.5) = 1.63134 M(-.25;1;-5) = 1.67206 M(-.25;1;-5.5) = 1.71009 M(-.25;1;-6) = 1.74579 M(-.25;1;-6.5) = 1.77947 M(-.25;1;-7) = 1.81136 M(-.25;1;-7.5) = 1.84167 M(-.25;1;-8) = 1.87056 M(-.25;1;-8.5) = 1.89818 M(-.25;1;-9) = 1.92466 M(-.25;1;-9.5) = 1.95009 M(-.25;1;-10) = 1.97456 Which would give table values of: static const float table[21] = { 0.82157f, 0.91549f, 0.99482f, 1.06298f, 1.12248f, 1.17515f, 1.22235f, 1.26511f, 1.30421f, 1.34025f, 1.37371f, 1.40495f, 1.43428f, 1.46195f, 1.48815f, 1.51305f, 1.53679f, 1.55948f, 1.58123f, 1.60212f, 1.62223f}; There is also a formula which asymptotically approaches M(a;b;-x) for high values of x that is: (gamma(b)/gamma(b-a))*(x^-a) Which for M(-.25;1;-x) is sqrt(sqrt(x))*1.10326 at x=10 this gives a value of 1.96191 (vs. what I think is the true value of 1.97456). The reason I've gotten into all this is I'm trying to vectorize with SSE intrinsics some of the slower loops in the preprocessor, and the hypergeom_gain function is the only thing stopping me from removing all the branches in the loops. I don't know how critical the accuracy of the function is to the performance of the preprocessor, but if that aforementioned approximation was good enough for all the values of x, it would really speed the loops up. Let me know what you think. I hope I'm helping out here (and not just confusing things). John Ridges Jean-Marc Valin wrote:> OK, so the problem is that the table you see does not match the definition? > y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x) > Note that the table data has an interval of .5 for the x axis. > > How far are your results from the data in the table? > > Cheers, > > Jean-Marc > > Quoting John Ridges <jridges at masque.com>: > > >> Thanks for the confirmation Jean-Marc. I kind of suspected from the >> comments that it was the confluent hypergoemetric function, which I was >> trying to evaluate using Kummer's equation, namely: >> >> M(a;b;x) is the sum from n=0 to infinity of (a)n*x^n / (b)n*n! >> where (a)n = a(a+1)(a+2) ... (a+n-1) >> >> But when I use Kummer's equation, I don't get the values in the >> "hypergeom_gain" table. Did you use a different solution to the >> confluent hypergoemetric function when you created the table? >> >> John Ridges >> >> >> Jean-Marc Valin wrote: >> >>> Hi John, >>> >>> M(;;) is the confluent hypergeometric function. >>> >>> Jean-Marc >>> >>> John Ridges a ?crit : >>> >>> >>>> Hi, >>>> >>>> I've been trying to re-create the table in the function "hypergeom_gain" >>>> in preprocess.c, and I just simply can't get the same values. I get the >>>> same value for the first element, so I know I'm computing gamma(1.25)^2 >>>> correctly, but I can't get the same numbers for M(-.25;1;-x), which I >>>> assume is Kummer's function. Is it possible that the comment is out of >>>> date and the values of Kummer's function used to make the table were >>>> different? Any help would be greatly appreciated. >>>> >>>> John Ridges >>>> >>>> >>>> _______________________________________________ >>>> Speex-dev mailing list >>>> Speex-dev at xiph.org >>>> http://lists.xiph.org/mailman/listinfo/speex-dev >>>> >>>> >>>> >>>> >>> >>> >> >> > > > > >
Jean-Marc Valin
2009-Jul-22 21:18 UTC
[Speex-dev] A technical question about the speex preprocessor.
Something looks odd without your values (or the doc) because hypergeom_gain() should really approach 1 as x goes to infinity. But in the end, an approximation is probably OK because denoising is anything but an exact science :-) Jean-Marc Quoting John Ridges <jridges at masque.com>:> By my reckoning the confluent hypergoemetric functions should have the > following values: > > M(-.25;1;-.5) = 1.11433 > M(-.25;1;-1) = 1.21088 > M(-.25;1;-1.5) = 1.29385 > M(-.25;1;-2) = 1.36627 > M(-.25;1;-2.5) = 1.43038 > M(-.25;1;-3) = 1.48784 > M(-.25;1;-3.5) = 1.53988 > M(-.25;1;-4) = 1.58747 > M(-.25;1;-4.5) = 1.63134 > M(-.25;1;-5) = 1.67206 > M(-.25;1;-5.5) = 1.71009 > M(-.25;1;-6) = 1.74579 > M(-.25;1;-6.5) = 1.77947 > M(-.25;1;-7) = 1.81136 > M(-.25;1;-7.5) = 1.84167 > M(-.25;1;-8) = 1.87056 > M(-.25;1;-8.5) = 1.89818 > M(-.25;1;-9) = 1.92466 > M(-.25;1;-9.5) = 1.95009 > M(-.25;1;-10) = 1.97456 > > Which would give table values of: > > static const float table[21] = { > 0.82157f, 0.91549f, 0.99482f, 1.06298f, 1.12248f, 1.17515f, 1.22235f, > 1.26511f, 1.30421f, 1.34025f, 1.37371f, 1.40495f, 1.43428f, 1.46195f, > 1.48815f, 1.51305f, 1.53679f, 1.55948f, 1.58123f, 1.60212f, 1.62223f}; > > There is also a formula which asymptotically approaches M(a;b;-x) for > high values of x that is: > > (gamma(b)/gamma(b-a))*(x^-a) > > Which for M(-.25;1;-x) is sqrt(sqrt(x))*1.10326 > > at x=10 this gives a value of 1.96191 (vs. what I think is the true > value of 1.97456). > > The reason I've gotten into all this is I'm trying to vectorize with SSE > intrinsics some of the slower loops in the preprocessor, and the > hypergeom_gain function is the only thing stopping me from removing all > the branches in the loops. I don't know how critical the accuracy of the > function is to the performance of the preprocessor, but if that > aforementioned approximation was good enough for all the values of x, it > would really speed the loops up. > > Let me know what you think. I hope I'm helping out here (and not just > confusing things). > > John Ridges > > > Jean-Marc Valin wrote: > > OK, so the problem is that the table you see does not match the definition? > > y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x) > > Note that the table data has an interval of .5 for the x axis. > > > > How far are your results from the data in the table? > > > > Cheers, > > > > Jean-Marc > > > > Quoting John Ridges <jridges at masque.com>: > > > > > >> Thanks for the confirmation Jean-Marc. I kind of suspected from the > >> comments that it was the confluent hypergoemetric function, which I was > >> trying to evaluate using Kummer's equation, namely: > >> > >> M(a;b;x) is the sum from n=0 to infinity of (a)n*x^n / (b)n*n! > >> where (a)n = a(a+1)(a+2) ... (a+n-1) > >> > >> But when I use Kummer's equation, I don't get the values in the > >> "hypergeom_gain" table. Did you use a different solution to the > >> confluent hypergoemetric function when you created the table? > >> > >> John Ridges > >> > >> > >> Jean-Marc Valin wrote: > >> > >>> Hi John, > >>> > >>> M(;;) is the confluent hypergeometric function. > >>> > >>> Jean-Marc > >>> > >>> John Ridges a ?crit : > >>> > >>> > >>>> Hi, > >>>> > >>>> I've been trying to re-create the table in the function "hypergeom_gain" > >>>> in preprocess.c, and I just simply can't get the same values. I get the > >>>> same value for the first element, so I know I'm computing gamma(1.25)^2 > >>>> correctly, but I can't get the same numbers for M(-.25;1;-x), which I > >>>> assume is Kummer's function. Is it possible that the comment is out of > >>>> date and the values of Kummer's function used to make the table were > >>>> different? Any help would be greatly appreciated. > >>>> > >>>> John Ridges > >>>> > >>>> > >>>> _______________________________________________ > >>>> Speex-dev mailing list > >>>> Speex-dev at xiph.org > >>>> http://lists.xiph.org/mailman/listinfo/speex-dev > >>>> > >>>> > >>>> > >>>> > >>> > >>> > >> > >> > > > > > > > > > > > > >
John Ridges
2009-Jul-22 22:18 UTC
[Speex-dev] A technical question about the speex preprocessor.
I got the approximation from a Google book: http://books.google.com/books?id=2CAqsF-RebgC&pg=PA385 Page 392, formula (10.33) Using this formula, you're right, hypergeom_gain() would *not* converge to 1 for large x, but would instead be gamma(1.25)/sqrt(sqrt(x)) which would approach zero. Now if the formula for the hypergeometric gain were instead gamma(1.5) * M(-.5;1;-x) / sqrt(x) that *would* approach 1, but that's just me noodling around with the formula to get something that approaches 1. Since I don't know how the hypergoemetric gain was derived (or even really what it means) I don't know if that's useful or not. Can you tell me what the source was for the original table values? John Ridges Jean-Marc Valin wrote:> Something looks odd without your values (or the doc) because hypergeom_gain() > should really approach 1 as x goes to infinity. But in the end, an > approximation is probably OK because denoising is anything but an exact science > :-) > > Jean-Marc > > Quoting John Ridges <jridges at masque.com>: > > >> By my reckoning the confluent hypergoemetric functions should have the >> following values: >> >> M(-.25;1;-.5) = 1.11433 >> M(-.25;1;-1) = 1.21088 >> M(-.25;1;-1.5) = 1.29385 >> M(-.25;1;-2) = 1.36627 >> M(-.25;1;-2.5) = 1.43038 >> M(-.25;1;-3) = 1.48784 >> M(-.25;1;-3.5) = 1.53988 >> M(-.25;1;-4) = 1.58747 >> M(-.25;1;-4.5) = 1.63134 >> M(-.25;1;-5) = 1.67206 >> M(-.25;1;-5.5) = 1.71009 >> M(-.25;1;-6) = 1.74579 >> M(-.25;1;-6.5) = 1.77947 >> M(-.25;1;-7) = 1.81136 >> M(-.25;1;-7.5) = 1.84167 >> M(-.25;1;-8) = 1.87056 >> M(-.25;1;-8.5) = 1.89818 >> M(-.25;1;-9) = 1.92466 >> M(-.25;1;-9.5) = 1.95009 >> M(-.25;1;-10) = 1.97456 >> >> Which would give table values of: >> >> static const float table[21] = { >> 0.82157f, 0.91549f, 0.99482f, 1.06298f, 1.12248f, 1.17515f, 1.22235f, >> 1.26511f, 1.30421f, 1.34025f, 1.37371f, 1.40495f, 1.43428f, 1.46195f, >> 1.48815f, 1.51305f, 1.53679f, 1.55948f, 1.58123f, 1.60212f, 1.62223f}; >> >> There is also a formula which asymptotically approaches M(a;b;-x) for >> high values of x that is: >> >> (gamma(b)/gamma(b-a))*(x^-a) >> >> Which for M(-.25;1;-x) is sqrt(sqrt(x))*1.10326 >> >> at x=10 this gives a value of 1.96191 (vs. what I think is the true >> value of 1.97456). >> >> The reason I've gotten into all this is I'm trying to vectorize with SSE >> intrinsics some of the slower loops in the preprocessor, and the >> hypergeom_gain function is the only thing stopping me from removing all >> the branches in the loops. I don't know how critical the accuracy of the >> function is to the performance of the preprocessor, but if that >> aforementioned approximation was good enough for all the values of x, it >> would really speed the loops up. >> >> Let me know what you think. I hope I'm helping out here (and not just >> confusing things). >> >> John Ridges >> >> >> Jean-Marc Valin wrote: >> >>> OK, so the problem is that the table you see does not match the definition? >>> y = gamma(1.25)^2 * M(-.25;1;-x) / sqrt(x) >>> Note that the table data has an interval of .5 for the x axis. >>> >>> How far are your results from the data in the table? >>> >>> Cheers, >>> >>> Jean-Marc >>> >>> Quoting John Ridges <jridges at masque.com>: >>> >>> >>> >>>> Thanks for the confirmation Jean-Marc. I kind of suspected from the >>>> comments that it was the confluent hypergoemetric function, which I was >>>> trying to evaluate using Kummer's equation, namely: >>>> >>>> M(a;b;x) is the sum from n=0 to infinity of (a)n*x^n / (b)n*n! >>>> where (a)n = a(a+1)(a+2) ... (a+n-1) >>>> >>>> But when I use Kummer's equation, I don't get the values in the >>>> "hypergeom_gain" table. Did you use a different solution to the >>>> confluent hypergoemetric function when you created the table? >>>> >>>> John Ridges >>>> >>>> >>>> Jean-Marc Valin wrote: >>>> >>>> >>>>> Hi John, >>>>> >>>>> M(;;) is the confluent hypergeometric function. >>>>> >>>>> Jean-Marc >>>>> >>>>> John Ridges a ?crit : >>>>> >>>>> >>>>> >>>>>> Hi, >>>>>> >>>>>> I've been trying to re-create the table in the function "hypergeom_gain" >>>>>> in preprocess.c, and I just simply can't get the same values. I get the >>>>>> same value for the first element, so I know I'm computing gamma(1.25)^2 >>>>>> correctly, but I can't get the same numbers for M(-.25;1;-x), which I >>>>>> assume is Kummer's function. Is it possible that the comment is out of >>>>>> date and the values of Kummer's function used to make the table were >>>>>> different? Any help would be greatly appreciated. >>>>>> >>>>>> John Ridges >>>>>> >>>>>> >>>>>> _______________________________________________ >>>>>> Speex-dev mailing list >>>>>> Speex-dev at xiph.org >>>>>> http://lists.xiph.org/mailman/listinfo/speex-dev >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>> >>>> >>> >>> >>> >>> >> >> > > > > >