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 >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>> >>>> >>> >>> >>> >>> >> >> > > > > >
Jean-Marc Valin
2009-Jul-23 00:07 UTC
[Speex-dev] A technical question about the speex preprocessor.
It's been a while since I did the maths. M(-.5;1;-x) optimises something else, though it's not too far. I think it may be [M(-.25;1;-x)]^2 (or is it the sqrt?) that is supposed to be there. In any case, if there's a mismatch between the doc and the code, the code is the one likely to be correct. Jean-Marc John Ridges a ?crit :> 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 >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>>>>> >>>>> >>>> >>>> >>>> >>>> >>> >>> >> >> >> >> >> > > >
John Ridges
2009-Jul-23 00:41 UTC
[Speex-dev] A technical question about the speex preprocessor.
That's got to be the answer. Using [M(-.25;1;-x)]^2 not only gives me exactly the values in the table, but it also makes gamma(1.25)^2 * [M(-.25;1;-x)]^2 / sqrt(x) converge to exactly 1 as x goes to infinity. Mystery solved. Thanks for all your help. John Ridges Jean-Marc Valin wrote:> It's been a while since I did the maths. M(-.5;1;-x) optimises something > else, though it's not too far. I think it may be [M(-.25;1;-x)]^2 (or is > it the sqrt?) that is supposed to be there. In any case, if there's a > mismatch between the doc and the code, the code is the one likely to be > correct. > > Jean-Marc > > John Ridges a ?crit : > >> 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 >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>>> >>>>>>> >>>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>>> >>>> >>>> >>> >>> >>> >>> >> >> > > >