[Speex-dev] A technical question about the speex preprocessor.
John Ridges
jridges at masque.com
Wed Jul 22 17:41:09 PDT 2009
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
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>
>>>
>>>
>>>
>>
>>
>
>
>
More information about the Speex-dev
mailing list