Fixed point arithmetic means assuming the binary number has an implicit decimal point at some place (for example, in the middle of a 32 bit value). Another way to think about this is to imagine that all stored numbers have been scaled by a factor (for example 2^16). One can even make it more complicated, and scale various parameters differently, in order to maximize the precision of the result for a given word size, but I won't discuss that here.
The beauty of fixed point arithmetic is that additions and subtractions are unaffected; one can use ordinary integer operations. Multiplications and divisions become more complicated, because they have to be scaled before or after the operation. Ideally one would make use of processor instructions such as double precision multiplies (that take two single precision numbers and produce a double precision result that can then be scaled). But to do this in C is tricky since you have to understand casting and promotion well, and have some idea about how to help the compiler optimization, since the language doesn't directly support mixed precision arithmetic.
For example, the following code does a multiply of two 32 values with an implicit 16 bit decimal point in C. Note the 16 bit shift right at the end to put the implicit decimal point back in the correct place. If compiled with optimization on for an x86 architecture, it will reduce down beautifully to a 32x32->64 bit multiply and a double-register shift.
Code:
#include <inttypes.h>
inline uint32_t mul(uint32_t i, uint32_t j)
{
return ((uint64_t) i* (uint64_t) j) >> 16;
}
Now to come back to your code; you are effectively trying to build an IIR filter. This is a simpler problem, particularly if you don't need the full 32 bit precision.
Basically you want to iteratively calculate
yi+1 = a.xi + b.yi. Typically in code one scales
y, so as to retain more precision; in your example, you scale
y by
10000 (which is called
cy), to give
cy = 10000*a*x + b*cy. The
a and
b are expressed as fractions that usually total to 1; the closer
a is to 1, the less effect the filter has.
What makes this an easier problem is that one can use the scaling to advantage so that one doesn't have to express the factors a and b as fractions. Instead of multiplying
x by
10000*a, instead use a new parameter
aa which is the fraction
a scaled by
10000, and treat it as an integer; now giving
cy = aa*x + b*cy. Similar, since
cy is scaled by
10000, then
b*cy is equivalent to
b*10000*cy/10000. So instead of multiplying
cy/10000 by
b*10000, instead use a new parameter
bb which is the fraction
b scaled by
10000, and treat it as an integer; now giving
cy = aa*x + bb*(cy/10000). The entire operation can now be done with integer arithmetic.
Ideally instead of 10000 one would choose a scaling factor that is a power of 2, so that the compiler can optimize the expensive divide into a shift. One can optimize it even further if
a is a binary fraction (removing all the multiplies), but I won't go into that here. One has to be careful that
cy does not overflow; this requires either the judicious use of double precision arithmetic, or simply use a scaling power that is half the length of the precision (eg, for 32 bits, use 65536 as the scaling factor, and keep all the parameters including
x within 16 bits, assuming unsigned arithmetic).
