Fixed Point Math

[poobah]

I'm just getting started in the XMOS world, and trying to implement a multi channel frequency generator using a XC-2 Ethernet kit. The goal is to generate varouis sine/cosine waves and mash them together in various combinations to produce spigrograph patterns for laser shows. I'm trying to implement this using fixed point math and sin lookup tables for the performance advantages. I'm using a 16:16 ratio that fits nicley in the 32-bit integer type.

My problem occurs when I try to multiply 2 numbers. I'm using the macs function to perform the multiply because it will likley overflow otherwise, but I'm stumped as how to shift the H and L registers back by 16 bits:

Code: Select all

int a=3;
int b=7;
int Fixed_a = a << 16;
int Fixed_b = b << 16;
{h,l} = macs(Fixed_a, Fixed_b, 0, 0);
Fixed_result = SomeMagicFunction (h,l);  
result = Fixed_result >> 16;  // should be 21 at this point

How should I shift H and L together in unison back by 16 bits in SomeMagicFunction? i.e. the LSB of H rolls into the MSB of L, and the LSB of L is discarded.
I'm sure that there's some asm code involved, but I haven't gotten that deep into the device yet.

A better question might be - has somebody done this already and is a library available? ;)

Thanks in advance,
Mark

[lilltroll]

If the values are limited to 16 bit's, using int32 and multiplication, it cannot overflow anyway.

What is you input range on a and b, [-2^15 2^15-1] ??
What 16:bit output do you want, the 16 MSB of that product ?

3*7 would end up in 0 in such case! (256*256 would be 1)

othervise

Code: Select all

{h, l} = macs(a, b , 0 , 0);
return ( (h<<1) + (l>>31) );

would return the 32 MSB of int32(a) * int32(b)

[poobah]

Thanks for the reply. I probably should have been a little more descriptive in my issue here.

Fixed point math shifts the decimal point over by a constant number of places, similar to an accountant's ledger. Instead of thinking of $137.48, you could represent this value as 13748 cents - Instead of 0.05, represent it as 5 percent, pi would be represented as 314, etc. This shift is accomplished by multiplying by a constant scalar value of 100.

Adding two of these number together is the same as dealing with decimal numbers. 0.05 + 1.32 = 1.37 and 5 + 132 = 137. (a*s) + (b*s) = (a+b)*s.

Multiplication becomes a problem because the scalar gets multiplied in twice - (a*s) * (b*s) = (a*b)*(s*s). In the example above, 13748 cents ($137.48) * 5 percent (0.05) becomes 68740 cents ($687.40) :shock: . The result needs to be divided by s again to become 687 cents ($6.87). the .4 of a cent can't be represented in the whole number and is lost.

Since computers deal with powers of 2 better than powers of 10, I'm doing a bitwise shift. In the 32 bit integer, I'm using a 16:16 ratio of whole number:fractional number (the scalar s is 65536). With this signed, it gives me a whole range of [-2^15 2^15-1] with an accuracy of 1/65536. The bits lay out like:

Code: Select all

Bit Number             32     31     30   ...   18     17   :    16     15   ...   2     1
Weight                Sign   2^15   2^14       2^1    2^0   :   2^-1   2^-2      2^-15 2^-16

2^-1 is 1/2
2^-2 is 1/4
2^-3 is 1/8 and so on

I could get greater resolution by changing the ratio to something like 12:20, but it would reduce the range of the whole portion.

pi would be represented as 0x0003:243F (the : is just to visually isolate the whole:fractional part).


So, I have the multiply working with my 32 bit numbers, but I need to shift the entire result back 16 bits. If I multiply pi above by 1.0 with the macs function it would look like (H and L in the product):

0003:243F * 0001:0000 = 0000:0003 2434:0000

If I just shift each back to the right independently by 16, the 3 rolls off the H value and L becomes 0000:2434.

In pure brute force C, i would do something like:

Code: Select all

int carry;
for (i=0;i<=16;i++)
{
  carry = 0;
  if (H && 0x00000001)
    carry = 1;
  H >> 1;
  L >> 1;
  if (carry==1)
    L = L || 0x80000000;
}
// check for overflow
if (H >0)
  L= 0xFFFFFFFF;

My final result would be in L. Anything that gets rolled off of that is smaller than 1/65535 and lost. Anything left in H represents an overflow. In assembly, I should be able to take advantage of a carry flag that makes this a lot more streamlined.

Sorry for the long post - i hope it answers your questions on what I'm trying to do here.

Thanks again!
Mark

[lilltroll]

OK

If a or b has at least 6 dB headroom, it's much more effective to use:

Assume b here,

Code: Select all

{h, l} = macs(a, b<<1 , 0 , 0);
return  h;

Would return the 32 MSB of int32(a) * int32(b)

The reason to the shift is just a consequence of the signed math, e.g.

(singed 31 bit ) * (singed 31 bit ) = ( signed 62 bit) == signed 30 bit + unsigned 32 bit
thus the 32 MSB is not bit 1-32 it's 2-33

[waluigi]

Code: Select all

{h, l} = macs(a, b<<1 , 0 , 0);
 :lol: return  h;

The reason to the shift is just a consequence of the signed math, e.g.

(singed 31 bit ) * (singed 31 bit ) = ( signed 62 bit) == signed 30 bit + unsigned 32 bit
thus the 32 MSB is not bit 1-32 it's 2-33

Nope, 1:0(16:16) * 1:0(16:16) = 0x10000 * 0x10000 = 0x100000000 = 1:0(32:32)

Code: Select all

int multi_1616 (int a, int b)
{
  int h, l;
  {h, l} = macs(a, b , 0 , 0);
  return  (h << 16) | (l >> 16);
}

[lilltroll]

[quote="waluigi"]

Code: Select all

int multi_1616 (int a, int b)
{
  int h, l;
  {h, l} = macs(a, b , 0 , 0);
  return  (h << 16) | (l >> 16);
}
[/quote]

Are you sure that you do not need 
unsigned l;  :?: 

otherwise the l>>16 will compile to ARSHI (arithmetic shift right immediate)  :!:  :?:  (I havn't tested your example in the compiler)

[lilltroll]

Seen from the computational burden, why use 16:16 ?

If you scale to {max min} |1 -1] , you get rotate 32, which doesn't need any shifting on a 32 bit machine, since it ends up in separate registers*.
I guess the amount of used instruction is of interest to you !?

The signed version might need the <<1

[lilltroll]

What about: :?:

-1:0(16:16) * 1:0(16:16) = 0xFFFF0000 * 0x00010000 = 0xFFFF00000000 = { you need rotate, not shift ?? }

-1:0(16:16) * -1:0(16:16) = 0xFFFF0000 * 0xFFFF0000 = 0xFFFE000100000000 = ???? != 1:0(32:32)

Is it:

Code: Select all

unsigned multi_1616 (unsigned a, unsigned b)
{
  unsigned h, l;
  {h, l} = mac(a, b , 0 , 0);
  return  (h << 16) | (l >> 16);
}

for the case with the frequency ?

[waluigi]

-1:0(16:16) * 1:0(16:16) = 0xFFFF0000 * 0x00010000 = 0xFFFF00000000 = { you need rotate, not shift ?? }

-1:0(16:16) * -1:0(16:16) = 0xFFFF0000 * 0xFFFF0000 = 0xFFFE000100000000 = ???? != 1:0(32:32)

In the first case, h = 0xFFFF and l = 0, so return 0xFFFF0000, correctly.
In the second case, h = 0xFFFE0001 and l = 0 so return 0x00010000, correctly.

BTW you might be right about the unsigned 'l'.

[lilltroll]

-1:0(16:16) * 1:0(16:16) = 0xFFFF0000 * 0x00010000 = 0xFFFF00000000 = { you need rotate, not shift ?? }

-1:0(16:16) * -1:0(16:16) = 0xFFFF0000 * 0xFFFF0000 = 0xFFFE000100000000 = ???? != 1:0(32:32)

In the first case, h = 0xFFFF and l = 0, so return 0xFFFF0000, correctly.
In the second case, h = 0xFFFE0001 and l = 0 so return 0x00010000, correctly.

BTW you might be right about the unsigned 'l'.

SHL Shift left
Shifts a word left by y bits, filling the least significant y bits with zeros.

-1:0(16:16) * -1:0(16:16) = 0xFFFF0000 * 0xFFFF0000 = 0xFFFE000100000000 =>

0xFFFF0001 <<16 0x0001 = 1

edited the post due to Copy&Paste error