0

I have found that piece of code on blackfin533, which has fract32 which is from -1, 1, its in the format 1.31.

I can't get why the pre-shifting is required for calculating the amplitude of a complex number (re, img). I know if you want to multiply 1.31 by 1.31 fractional format then you need to shift right 31 bits.

GO_coil_D[0].re, and GO_coil_D[0].im are two fract32.

I can't get what the following code is doing :

 norm[0] = norm_fr1x32(GO_coil_D[0].re);
 norm[1] = norm_fr1x32(GO_coil_D[0].im);
 shift = (norm[0] < norm[1]) ? (norm[0] - 1) : (norm[1] - 1);
 vectorFundamentalStored.im = shl_fr1x32(GO_coil_D[0].im,shift);     
 vectorFundamentalStored.re = shl_fr1x32(GO_coil_D[0].re,shift);
 vectorFundamentalStored.im = mult_fr1x32x32(vectorFundamentalStored.im, vectorFundamentalStored.im);
 vectorFundamentalStored.re = mult_fr1x32x32(vectorFundamentalStored.re, vectorFundamentalStored.re);  
 amplitudeFundamentalStored = sqrt_fr16(round_fr1x32(add_fr1x32(vectorFundamentalStored.re,vectorFundamentalStored.im))) << 16;
 amplitudeFundamentalStored = shr_fr1x32(amplitudeFundamentalStored,shift);

round_fr1x32` (fract32 f1) fract16 Rounds the 32-bit fract to a 16-bit fract using biased rounding.

norm_fr1x32 norm_fr1x32 (fract32) int Returns the number of left shifts required to normalize the input variable so that it is either in the interval 0x40000000 to 0x7fffffff, or in the interval 0x80000000 to 0xc0000000. In other words, fract32 x; shl_fr1x32(x,norm_fr1x32(x)); returns a value in the range 0x40000000 to 0x7fffffff, or in the range 0x80000000 to 0xc0000000

andre_lamothe
  • 2,171
  • 2
  • 41
  • 74

1 Answers1

1

1) If the most significant n bits of the fractional part are all '0' bits, and they are followed by a '1' bit, then n behaves like a floating point binary exponent of value n, and the remaining 31-n bits behave like the mantissa. Squaring the number doubles the number of leading '0' bits to 2*n and reduces the size of the mantissa to 31-2*n bits. This can lead to a loss of precision in the result of the squaring operation.

2) round_fr1x32 converts the 1.31 fraction to a 1.15 fraction, losing up to 16 more bits of precision.

Hopefully you can see that steps 1 and 2 can remove a lot of precision in the number. Pre-scaling the number reduces the number of leading '0' bits n as much as possible, resulting in less precision being lost at step 1. In fact, for one of the two numbers being squared and added, the number of leading '0' bits n will be zero, so squaring that number will still leave up to 31 bits of precision before it is added to the other number. (Step 2 will reduce that precision to 15 bits.)

Lastly, you are wrong about the result of multiplying two 1.31 fraction format numbers - the result needs to be shifted right by 31 bits, not 62 bits.

Worked example:

Let's say the real part is 3/1024 and the imaginary part is 4/1024 in decimal, so the absolute value should be 5/1024 by pythagoras.

With no pre-scaling, the binary fractions are re=0.0000000011₂, im=0.0000000100₂. Squaring them gives re²=0.00000000000000001001₂, im²=0.00000000000000010000₂. Adding the squares gives abs²=0.00000000000000011001₂. Rounding to 15 fractional bits gives abs²=0.000000000000001₂. Taking the square root gives abs=0.000000010110101₂. This differs from the exact result 0.0000000101₂ by 0.000000000010101₂.

When prescaling, both fractions are shifted left by 6 bits, giving sre=0.0011₂, sim=0.0100₂ (I used the prefix 's' to mean 'scaled'). Squaring them gives sre²=0.00001001₂, sim²=0.00010000₂. Adding the squares gives sabs²=0.00011001₂. Rounding to 15 fractional bits does not change the value. Taking the square root gives sabs=0.01010000₂. Converting that to 1.31 format and shifting right by 6 bits gives abs=0.0000000101₂ which is exactly correct (5/1024 in decimal).

Ian Abbott
  • 15,083
  • 19
  • 33
  • 1
    I added a worked example with and without prescaling for the case when the real part is 3/1024 (in decimal) and the imaginary part is 4/1024. Hopefully you can see that the absolute value (amplitude) should be 5/1024 since it forms a pythagoras 3-4-5 triangle. – Ian Abbott Oct 03 '16 at 17:51
  • Would you tell me what's the original problem of multiplying the real with itself and the imaginary with itself, then adding the result, then take the square root, then shift right 31 bits ? what would be the problem with that approach ? – andre_lamothe Oct 03 '16 at 18:49
  • 1
    If you can calculate the square root of a 1.31 format fraction easily, there is not much problem doing it that way, although there is likely to be some loss of bits of precision when the real and imaginary parts are squared. – Ian Abbott Oct 04 '16 at 09:32
  • how would I calculate the square root of a 1.31 format ? – andre_lamothe Oct 04 '16 at 10:11
  • there is a function in the library that calculates the square root of fract32, is that the one that should I use ? or should I do some pre conditioning of the inputs ? – andre_lamothe Oct 04 '16 at 10:17
  • 1
    There is no harm in trying the sqrt_fr32 function. Then the line `amplitudeFundamentalStored = sqrt_fr16(round_fr1x32(add_fr1x32(vectorFundamentalStored.re,vectorFundamentalStored.im))) << 16;` would change to `amplitudeFundamentalStored = sqrt_fr32(add_fr1x32(vectorFundamentalStored.re,vectorFundamentalStored.im));`. You might still want to do the shifting for normalization if you want the result to be as precise as possible. – Ian Abbott Oct 04 '16 at 13:05
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/124898/discussion-between-ahmed-saleh-and-ian-abbott). – andre_lamothe Oct 04 '16 at 13:09