fixed-point approximation

Can anyone come up with a 'good' fixed-point approximation for the
follow mapping

Q = sqrt(2^B) / (2^B - 1)

where B is in range 0.1 to 10. Only integer arithmetic is allowed. B
could for instance be represented in Q10.

The expression maps a bandwidth measured in octaves to a Q factor
(center frequency divided by bandwidth).

On 6/27/2011 8:28 AM, niarn wrote:
> Can anyone come up with a 'good' fixed-point approximation for the
> follow mapping
>
> Q = sqrt(2^B) / (2^B - 1)
>
> where B is in range 0.1 to 10. Only integer arithmetic is allowed. B
> could for instance be represented in Q10.
>
> The expression maps a bandwidth measured in octaves to a Q factor
> (center frequency divided by bandwidth).
>

Couple things that make this question odd:  'Q' is overloaded.  You are 
using it to specify the result of the equation and fixed-point type and 
Q10 is ambiguous.  I assume Q10 you mean 10 fraction bits and one sign bit?

 From my perspective, define the *range* and *resolution* for B and you 
have your requirements.  You defined your range 0.1 - 10, then determine 
the resolution (smallest fractional required) then you can determine 
what fixed point representation needed, the number of fractional bits 
and number of integer bits.

Chris

niarn wrote:

> Can anyone come up with a 'good' fixed-point approximation for the
> follow mapping
> 
> Q = sqrt(2^B) / (2^B - 1)
> 
> where B is in range 0.1 to 10. Only integer arithmetic is allowed. B
> could for instance be represented in Q10.
> 
> The expression maps a bandwidth measured in octaves to a Q factor
> (center frequency divided by bandwidth).

What you need is a 2^x approximation, where x = [0...1]. Everything else 
is trivial arithmetic. Depending on the accuracy requirements, the 2^x 
could probably be done as simple as a 2nd order approximation.


Vladimir Vassilevsky
DSP and Mixed Signal Design Consultant
http://www.abvolt.com

Sorry, I just realized that there are answers to my post. I have been using
google groups until now but there I can't see any activity.

@Christopher Felton
Sorry, for the confusion. I was hoping that the meaning of Q and Q10 would
be clear from the context, but obviously not. Yes, Q10 means 10 bits for
the fractional part. I don't care if there is a sign bit but there doesn't
have to be one.

The hard part is not in determining an appropriate representation for the
input (B) and output (Q). The mapping describes a monotonically decreasing
function, for B=[0.1, 10], it has a shape similar to (1/x). I tried
experimentally to approximate the mapping using 1/x and exp(-x) but could
not come up with something good. I wasn't satisfied with a 5th order Taylor
expansion around B=1. So after 1 to 2 hours of fun I finally made a LUT.
The reason why I posted the 'challenge' here is that the mapping has some
relevance to (practical) DSP and it is maybe not so straight forward to
find a 'good' and 'cheap' approximation if LUT approach is not considered.

@Vladimir Vassilevsky
I did not consider LUT for 2^x because it's more direct to make LUT for the
mapping itself.

Cheers

looks like google groups is working again.

On Jun 27, 9:28&#4294967295;am, niarn <niar...@gmail.com> wrote:
> Can anyone come up with a 'good' fixed-point approximation for the
> follow mapping
>
> Q = sqrt(2^B) / (2^B - 1)
>

looks like you got this out of the Wikipedia page on equalization
filters.  (it relates Q and bandwidth in octaves.)

it's equivalent to the equation in the Audio EQ Cookbook, but does not
have the compensation of the BLT frequency warping (that scrunches
bandwidth).

> where B is in range 0.1 to 10. Only integer arithmetic is allowed. B
> could for instance be represented in Q10.
>
> The expression maps a bandwidth measured in octaves to a Q factor
> (center frequency divided by bandwidth).

if you have division available, it's a simple sinh() function.

   1/Q = 2*sinh( ln(2)/2 * B )

or

   Q = 0.5/sinh( ln(2)/2 * B )



this has a power series with odd-order terms that converge quickly.
just pick out the first 2 or 3 terms.

or the inverted function can have coefficients for the odd-order terms
solved directly with the Remez exchange alg.  i might be able to run
that in Octave.  i can send you code if you want.

r b-j

On Jun 27, 9:28&#4294967295;am, niarn <niar...@gmail.com> wrote:
> Can anyone come up with a 'good' fixed-point approximation for the
> follow mapping
>
> Q = sqrt(2^B) / (2^B - 1)
>
> where B is in range 0.1 to 10. Only integer arithmetic is allowed. B
> could for instance be represented in Q10.
>
> The expression maps a bandwidth measured in octaves to a Q factor
> (center frequency divided by bandwidth).

Try 1.44269504/B  - 0.02888113B

You may round the coefs to get your fixed point form, plus the
reciprocal may be found iteratively quite efficiently - think Newton-
Raphson. There are other ways too.

IHTH,
Clay

Thanks r b-j. Great response.

An approximation based on a taylor expansion of sinh gives good results.
The following approximation

x = log(2)/2 * B;
y = x + (1/6)*x.^3 + (1/120)*x.^5 + (1/5040)*x.^7;
Q = 0.5./y;

gives a maximum relative error of ~0.05% at B = 6, when the interval
B=[0.2, 6] is considered.

The following implementation gives a maximum relative error of ~0.1% and
works for B=[0.2, 6].
Input is in Q10 and output is returned in Q12.

#define LOG2DIV2Q15    11357
#define ONEOVER6Q33    1431655765    
#define ONEOVER120Q37  1145324612    
#define ONEOVER5040Q43 1745256552

short B2Q(short sB10)
{
	long long x, x2, x3, x5, x7;
	int acc;
	long long lTmp;

	x  = LOG2DIV2Q15*sB10; 
	acc = (int)(x << 4);   
	lTmp = x*x;			
	x2 = (lTmp >> 22);   
	lTmp = x*x2;
	x3 = (int)(lTmp >> 26); 
	lTmp = x3*ONEOVER6Q33;
	acc += (int)(lTmp >> 31); 
	lTmp = x2*x3;
	x5  =  lTmp >> 30; 
	lTmp = x5*ONEOVER120Q37;
	acc += (int)(lTmp >> 33); 
	lTmp = x5*x2;
	x7 = lTmp >> 30; 
	lTmp = ONEOVER5040Q43*x7;
	acc += (int)(lTmp >> 37);  
	lTmp = (1LL << 56) / acc;  
	return (short)(lTmp >> 16);
}

r b-j, forgot to mention that I got the expression for the B to Q mapping
from here
http://www.astralsound.com/parametric_eq.htm
and not from wiki.

Btw. where did you get the expression involving sinh from? 

Thanks Clay. That is a very neat approximation. Initially I looked at an
approximation similar to this one. However, it's relative error increases
with B and at B=6 the relative error is ~50%. However if speed is crucial
(such as if the computation resides in a process loop) this approximation
can be handy. 

On Jul 1, 7:18=A0am, "niarn" <niaren9@n_o_s_p_a_m.gmail.com> wrote:
> r b-j, forgot to mention that I got the expression for the B to Q mapping
> from herehttp://www.astralsound.com/parametric_eq.htm
> and not from wiki.
>
> Btw. where did you get the expression involving sinh from?

i derived it a couple of decades ago.  first for analog BPF:

   1/Q =3D 2*sinh( ln(2)/2 * BW )

and then for the same filter mapped to digital using the bilinear
transform:

   1/Q =3D 2*sinh( ln(2)/2 * BW * w0/sin(w0) )

where w0 =3D 2*pi*f0/Fs is the center frequency (angular and normalized)

the analog form is equivalent to what you got from wherever site (and
is the same as in the wikipedia article).  the digital form has the BW
stretched out slightly by the 1/sinc factor to make a first-order
approximation to compensate for the scrunching of the bandwidth done
by the bilinear transform.

r b-j