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
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).