-
Given two real (floating point or fixedpoint), `a` and `b`, where `b` can be negative, how to efficiently compute the power function `pow(a, b)` ?
Mathematically, it's required to compute pow(a, b)=exp(b*ln(a))
I have to implement it in C, and then, I have to rewrite it in assembly :)
edit: Typo: a^b.
a,b are real numbers
b can be negative.
-
Aren't Taylor series and Newton-Raphson's method the answer to most computations on irrational functions?
Also, google is your friend
https://stackoverflow.com/questions/4518011/algorithm-for-powfloat-float
https://stackoverflow.com/questions/9799041/efficient-implementation-of-natural-logarithm-ln-and-exponentiation/44232045
-
Other than for a learning exercise, the answer normally is to use the compiler's floating/fixed point library and a powerful enough CPU, with sufficient memory, to handle it.
-
Given two real (floating point or fixedpoint), `a` and `b`, where `n` can be negative,
And 'n' would be? Is it 'a' or 'b'?
As it makes quite a complex difference (pun totally intended).
-
b being negative isn't an issue, that's just fractional division rather than multiplication, and the definition of the exponent function provides continuity. Surely, a being negative is the bigger problem -- fractional powers are imaginary but integer powers are just fine?
https://en.cppreference.com/w/c/numeric/math/pow
Tim
-
Given two real (floating point or fixedpoint), `a` and `b`, where `b` can be negative, how to efficiently compute the power function `pow(a, b)` ?
Mathematically, it's required to compute pow(a, b)=exp(b*ln(a))
I have to implement it in C :)
I suppose you're not allowed to use the C std math library, because it's all there already. And if you have access to a compliant C compiler, there's no reason you wouldn't have access to this library.
So I'm assuming it's just an exercise. Or maybe you're looking for something more efficient (in which case, have you benchmarked the standard pow() function already to know what you'd be competing against?)
You can take a look at the following implementation in openlibm:
https://github.com/JuliaMath/openlibm/blob/master/src/e_pow.c
-
You can take a look at the following implementation in openlibm:
https://github.com/JuliaMath/openlibm/blob/master/src/e_pow.c
That code doesn't look like a homework assignment. I can imagine a team of Applied Math PhDs sitting around discussing the boundary conditions for weeks. Then they spend even more time coming up with the algorithm for the programmer to code. There's some magic math in that code.
-
You can take a look at the following implementation in openlibm:
https://github.com/JuliaMath/openlibm/blob/master/src/e_pow.c
That code doesn't look like a homework assignment. I can imagine a team of Applied Math PhDs sitting around discussing the boundary conditions for weeks. Then they spend even more time coming up with the algorithm for the programmer to code. There's some magic math in that code.
Definitely. But it gives the basic calculations they use, so that can serve as an inspiration for implementing the same, even if you end up with something a lot less optimized.
Now I don't know if the OP's question was actually a homework assignment. If it was, then giving them a "homework assignment" directly usable answer would not do them any favor IMHO. And I also don't know what the value of such an assignment would be. If it were a programming assignment, then I don't really see the point of using anything else than the std pow() function. But if it were an applied math assignment, then it would be interesting, but in this case I would expect the students to at least have a rough idea of what kind of approach they can use to get an approximation of this function.
-
BKM might be useful. It is a bit difficult to grok though.
https://en.m.wikipedia.org/wiki/BKM_algorithm
-
just an exercise
Kind of personal challenge for which I will have to rewrite it in pure assembly in order to fit a very old CPU with limited resources.
Crazy, I know ;D
It's done for a larger hobby project, where for historical reasons I want to support all the technology I encountered during the last 25 years of my job and hobby experiences, so in the main big project there are several and different pieces of technology, including things used in the 80s, 90s, 2000s, 2010s and 2020s.
For the piece of the puzzle belonging to the 80s, I chose to use a very very old MPU to realize a height sensor, which is used in a part of the big project; in its planned scheme done on Labview as "proof of concept" between the CPU and the analogical sensor, Data linearization is done using the nonlinear curve fitting method, which requires three pow(a,b) calculations like:
h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q
a1,a2,a3 come from ADC sampling, so the final functionality will be height = function(ADC sampled value).
b1,b2,b3,k1,k2,k3 come from a wizard functionality offered by Labview to synthesize an interpolator.
The CPU is clocked at 2Mhz, it's very slow and limited compared to a modern MPU, there are no time-critical constraints but rather ram-rom-space constraints, since there is only 768 byte of RAM, 1Kbyte of UV-EPROM for the whole program, which should also do other stuff, like talking on the serial line, sampling data on two analog channels, properly filter data, and handling two digital lines to control a couple of motors.
The current firmware is already written in assembly. The pow(a,b) will be written in C, then translated into assembly, deep studied, and finally rewritten from scratch in assembly.
-
There is certainly a floating point library in asm (or binary) written for a similar MCU (ie from the same family but with larger rom) available.. For example a 32bit floating point lib (incl. sin, cos, tan, exp, ln, pow..) is usually ~3.5kB large with those 8bitters MCUs.. Even I doubt you would be able to cut a bigger chunk off the lib such it fits into your 1kB incl your other stuff..
In past 1kB was enough for a complete 8-9digits floating point math in vintage scientific calculators - but those MCUs were of a different arch (BCD instead of BIN)..
-
Does LabVIEW not produce integer coefficients (a simple polynomial fit)?
If not, you will find precisely because of this, why it is cheaper to do so. On a machine with hardware multiply, it's just a matter of pulling down parameters and chugging, on loop (Horner's method); without, it's calls/jumps to a multiplication routine (typically shift and accumulate). If you can constrain LabVIEW to produce results of this form, you will have a much easier -- and faster -- time.
There are also rational approximation methods ( https://en.wikipedia.org/wiki/Pad%C3%A9_approximant ), which on an add-shift machine, might be competitive with polynomial methods -- division is a shift-subtract algorithm, taking similar time to multiplication but fewer terms may be necessary.
I suppose this isn't an excuse to wholly avoid implementing pow() -- but it can push off the need for it quite far, and indeed the same method, loaded with different coefficients, might serve to implement that very function, at least over a certain range, and the range (and selection of coefficients) can then be managed by an outside function.
How many things are you looking to get into that MPU at once, anyway? And how much storage does it have accessible? If practically unlimited (e.g. parallel or SPI Flash chip, SD card*), how will you access and manage it? With so little ROM to work with, it seems likely that some kind of very basic, very general, interpreted language will be necessary. Perhaps even just to implement access to such things. Which means any real "thinking" (like evaluating pow()) will take very long indeed -- fast enough for a handheld calculator perhaps, but doubtful for something like a control loop. Well, maybe like temperature control, there's a lot that can be done between seconds of samples. But besides that...
*Probably not, at least not in the most convenient form -- FAT16 takes a few kB to implement I think, and is only good to a few hundred MB capacity (or a couple GB with some inconvenience); FAT32 takes a few kB more, and is good for pretty much whatever size is available these days. You might be able to access the SD card as a flat memory store, with less code, but not with FAT, and so you'll have to copy data to the SD card using some raw method (disk image, hex editor, etc.).
Hmm, I wonder if you can format an SD card for FAT32 and place exactly one file on it, and have that file be at a consistent sector offset, no fragmentation, etc...
Tim
-
Given the very limited resources, doing any FP might be extremely challenging - probably impossible. I'm assuming your old CPU doesn't have any FPU. Doing FP in software with this amount of memory? I doubt it.
Anyway, it's probably what T3sl4co1l more or less suggested, but I would switch to piecewise linear interpolation instead of trying to compute those power functions. Just find a decent piecewise linear approximation of each power term. Then the only computation needed will be segment lookup, linear interpolation and additions. All this using integer/fixed-point arithmetic.
-
What precision is your fixed or floating point? To what precision do you need the result? Do you have an FPU? How much program code size is OK to use? How much ROM is ok to use for tables? How much working RAM can you use?
For example, if you are using any 8 bit fixed point or floating point format, and you can afford 64 KB of ROM, then you can get the most accurate possible result with one shift-and-OR and a memory load. BOOM!
-
Does LabVIEW not produce integer coefficients (a simple polynomial fit)?
Unfortunately no. I get irrational coefficients.
-
there is only 768 byte of RAM, 1Kbyte of UV-EPROM for the whole program, which should also do other stuff
Ugh.
You're clearly not going to be aiming for a fast implementation. Just getting something to fit into that size at *all* will be challenging.
It would have been useful to have this information at the outset.
We still don't know what precision and range you need to work with, either for the inputs or the result.
-
To what precision do you need the result?
I am trying to amplify the signal 4X in order to achieve better precision. The on board ADC is an 8bit unit working in the range of 5V (VRefL= 0V VRef_H= 5V), which means ~20mV
5V/(255) is ~20mV
But! Since the measure range is 600..1300mm, I can isolate the corresponding sub voltage range, which is around 0V..1.2V, and can be 4X amplified by an OA.
Still ~20mV, but if the analogical sensor is stable enough, this will produce better numerical results.
Do you have an FPU?
No :D
How much program code size is OK to use?
At the moment I am using 768 byte of ram in total, but I am ready to prepare a voltage pump in order to program the onboard UV-Eprom (1Kbyte).
For example, if you are using any 8 bit fixed point or floating point format, and you can afford 64 KB of ROM, then you can get the most accurate possible result with one shift-and-OR and a memory load. BOOM!
64Kbyte is impossible for my system. The best I can do is 32Kbyte for both RAM and ROM.
-
To what precision do you need the result?
I am trying to amplify the signal 4X in order to achieve better precision. The on board ADC is an 8bit unit working in the range of 5V (VRefL= 0V VRef_H= 5V), which means ~20mV
5V/(255) is ~20mV
But! Since the measure range is 600..1300mm, I can isolate the corresponding sub voltage range, which is around 0V..1.2V, and can be 4X amplified by an OA.
Still ~20mV, but if the analogical sensor is stable enough, this will produce better numerical results.
But what numeric range do you consider these 8 bit values to be in? Is it the same for A and B? And the result?
Is B an arbitrary number, or is there only a limited set of possibilities that depend on the algorithm, not on the ADC data?
Depending on your answers, you can very quickly find that most of the results are either 0 or infinity.
-
But what numeric range do you consider these 8 bit values to be in?
the built-in ADC is 8bit, with an equivalent bit noise of 1 bit after filtering, that means what you read and filter is truly 20mV of resolution. Filtering costs something like 200 cycles due to filtering window.
The plan is converting the raw data from the ADC into floating-point/fixed-point, applies some interpolating function
Is it the same for A and B? And the result?
Values "a" and "b" in the exp(a,b) function come from Labview, they are both real rumbers, and they require to be floatingpoint or fixedpoint, but I cannot yet qualify how many fractional bits they need because I am still playing with Labview in order to optimize those values.
Is B an arbitrary number, or is there only a limited set of possibilities that depend on the algorithm, not on the ADC data?
There will be three "b" values, "b1", "b2", "b3", to be used in three "exp" function calls, but they are guaranteed to be constant (Labview says that), while "a" will be a kind of constant that multiply the value from the adc, kind of "a = k * adc_value"
and you have to calculate exp(a,b), so exp ( k_const * adc_value, b_const)
Depending on your answers, you can very quickly find that most of the results are either 0 or infinity.
Yup, including the LUT strategy, which seems more feasible although I haven't yet understood how to do it correctly ;D
-
This stuff is so impossible that I really wonder how Nasa did program the MPUs for the Apollo space missions. Their guidance computer (https://en.wikipedia.org/wiki/Apollo_Guidance_Computer) has 2K words of erasable magnetic-core memory and 36 kilowords of read-only core rope memory :o :o :o
And this stuff has to do more tasks than just controlling the position of a vertically movable platform ... always impressed by those Big Guys!
-
This stuff is so impossible that I really wonder how Nasa did program the MPUs for the Apollo space missions. Their guidance computer (https://en.wikipedia.org/wiki/Apollo_Guidance_Computer) has 2K words of erasable magnetic-core memory and 36 kilowords of read-only core rope memory :o :o :o
And this stuff has to do more tasks than just controlling the position of a vertically movable platform ... always impressed by those Big Guys!
I'm sure you must be overthinking all of this. It must be able to be linearised or something over the values in the possible range. And especially only talking about 8 bit precision anyway.
-
I'm sure you must be overthinking all of this. It must be able to be linearised or something over the values in the possible range. And especially only talking about 8 bit precision anyway.
How can you be sure if you haven't seen anything?
I wrote the above showed piece of code to acquire data and send them over the serial, so I sampled data within the range {15..1600}mm and put these samples inside a Labview's table, and when I click on the button it suggested a linearization model made with three pow(a,b), as I described above.
It means, if you want to linerize these
...
2.75 015 cm
2.55 020 cm
2.00 030 cm
1.55 040 cm
1.25 050 cm
1.15 060 cm
0.90 070 cm
0.80 080 cm
0.75 090 cm
0.65 100 cm
0.60 110 cm
0.55 120 cm
0.50 130 cm
0.44 140 cm
0.45 150 cm
...
You need to implement this: h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q
You have to optimize a1,b1,a2,b2,a3,b3,q, but you will always have to deal with fpow, fmul, fadd, fand fcmp
So I tried to consider only the subrange of my interest, just to see if it somehow look "locally linear".
1.15 060 cm
0.90 070 cm
0.80 080 cm
0.75 090 cm
0.65 100 cm
0.60 110 cm
0.55 120 cm
0.50 130 cm
Unfortunately, Labview says that it does not, and still suggest to synthesize two exp(a,b) functions.
So we have just moved from this
h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q
to this
h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + q
-
What is the precision/accuracy of the sensor?
How did you get the data?
This part of the data suggests there the sensor measurements are not "very accurate":
...
0.50 130 cm
0.45 150 cm
0.44 140 cm
...
My suggestion is transform x into 1/x:
0,36364 15
0,39216 20
0,5 30
0,64516 40
0,8 50
0,86957 60
1,11111 70
1,25 80
1,33333 90
1,53846 100
1,66667 110
1,81818 120
2 130
2,27273 140
This way you can get a more linear relation.
-
It means, if you want to linerize these
...
2.75 015 cm
2.55 020 cm
2.00 030 cm
1.55 040 cm
1.25 050 cm
1.15 060 cm
0.90 070 cm
0.80 080 cm
0.75 090 cm
0.65 100 cm
0.60 110 cm
0.55 120 cm
0.50 130 cm
0.44 140 cm
0.45 150 cm
...
You need to implement this: h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q
Maybe I'm misunderstanding something, but have you tried piecewise curve fitting?
For the above numbers you get a pretty good fit with quadratic curve fitting if you break up the input range into three intervals: [2.75 .. 1.55], [ 1.55 .. 1.15 ] and [ 1.15 .. 0.45 ].
-
Unfortunately the sensor doesn't have any good documentation.
I don't know anything about its precision/accuracy.
_______
/0..20Hz|
sensor _____/ OA G=1 |_____ ADC 8bit _____ raw data over serial line
\follower| CPU builtin
\_______|
This is what I am doing at the moment, and how I collected data to prepare the above tables.
but have you tried piecewise curve fitting?
Not yet. Good point :D
I have only tried the "NI myDAQ" curve fitting software offered by LabVIEW.
The CPU reads the ADC, and sends raw data to the computer without any manipulation.
The PC collects data and passes them to the myDAQ software.
I am also preparing a similar working scenario for Octave + GNUplot.
For the above numbers you get a pretty good fit with quadratic curve fitting if you break up the input range into three intervals: [2.75 .. 1.55], [ 1.55 .. 1.15 ] and [ 1.15 .. 0.45 ].
Nice! Thanks! Maybe breaking up the input range into several intervals is the best way to sort it out!
Thanks, again, I hadn't thought about it yet :D
-
I think I would look at a polynomial curve fit and graph the actual and interpolated data. I don't know anything about Labview but Octave (or MATLAB) can do this with a single statement to any degree required.
The link covers single polynomial fitting and spline fitting:
https://octave.org/doc/v4.4.1/Polynomial-Interpolation.html
If the computation is to be done on the uC, all I would need to store in memory are the coefficients.
Clearly, sending the raw data to the PC and using more powerful math tools will be a lot easier.
-
Just to consider alternatives, which sensor would you use to (contactless) measure a distance in the range of 60cm .. 130cm?
-
Just to consider alternatives, which sensor would you use to (contactless) measure a distance in the range of 60cm .. 130cm?
What kind of sensor are you considering in this thread? I may have missed it, but all I saw was "analog sensor". Just needs a bit of details so we don't suggest the same.
For this range, ultrasonic sensors may be the cheapest and easiest approach? 5mm accuracy should be doable.
Otherwise, you may consider time-of-flight sensors such as this: https://www.st.com/en/imaging-and-photonics-solutions/vl53l0x.html (https://www.st.com/en/imaging-and-photonics-solutions/vl53l0x.html) (not sure you'd manage 5mm accuracy with this though...)
-
What kind of sensor are you considering in this thread?
A sensor that without any contact can measure a distance in the range of 60cm .. 130cm.
I am currently using a pure analogical optical reflective IR sensor.
I really like the vl53l0x (https://www.st.com/en/imaging-and-photonics-solutions/vl53l0x.html) sensor. It's embedded, it's powerful, and it's very small. I think I will use one of these for other projects, thanks :D
-
I really like the vl53l0x (https://www.st.com/en/imaging-and-photonics-solutions/vl53l0x.html) sensor. It's embedded, it's powerful, and it's very small. I think I will use one of these for other projects,
Those are nice and reasonably cheap, but I'm not too sure they would yield the kind of accuracy you're looking for. Something to prototype though.
-
I'm sure you must be overthinking all of this. It must be able to be linearised or something over the values in the possible range. And especially only talking about 8 bit precision anyway.
How can you be sure if you haven't seen anything?
Experience.
I wrote the above showed piece of code to acquire data and send them over the serial, so I sampled data within the range {15..1600}mm and put these samples inside a Labview's table, and when I click on the button it suggested a linearization model made with three pow(a,b), as I described above.
It means, if you want to linerize these
...
2.75 015 cm
2.55 020 cm
2.00 030 cm
1.55 040 cm
1.25 050 cm
1.15 060 cm
0.90 070 cm
0.80 080 cm
0.75 090 cm
0.65 100 cm
0.60 110 cm
0.55 120 cm
0.50 130 cm
0.44 140 cm
0.45 150 cm
...
You need to implement this: h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q
What are the labels for the above columns? What is the input? What is the output? How does this relate to the formula? What on Earth are pow1, pow2, pow3?
Clearly your data are not very accurate given the reversal of 0.44 and 0.45.
You have to optimize a1,b1,a2,b2,a3,b3,q, but you will always have to deal with fpow, fmul, fadd, fand fcmp
If you are dealing with fpow at runtime on the microcontroller then this is not anything I know of as linearisation.
-
So we have just moved from this
h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + K3*pow3(a3,b3) + q
to this
h = K1*pow1(a1,b1) + K2*pow2(a2,b2) + q
Do you have values for K1 & K2, and the ranges and precision required for a1, b1, a2 and b2? I would like to play around with them...
-
What are the labels for the above columns? What is the input? What is the output?
The sensor outputs a voltage value for each distance between the sensor and the ground I try to measure. The measuring plane moves vertically, and I have a meter tape to take the measures. So I can report { Voltage, distance } in a table.
-
Have you considered an optical encoder -- either linear or rotary?
-
Do you have values for K1 & K2, and the ranges and precision required for a1, b1, a2 and b2? I would like to play around with them...
Abandoned for now since Labview keeps suggesting irrational numbers for those constants. I tried several times and always got the same result, an this seems to force me to use floating point in double precision, which is really too much for this project, and it's not worth it.
I am going to play with a quadratic approach. It seems more appropriate for my little CPU. Probably the National Instruments curve fitting assistant is made for a PC, where you have floating point in double precision for free as implemented in hardware and with large available resources (lot of ram, high clock frequency, etc).
None of this makes any sense.
Coefficients spat out by a program can be implemented by floating point or fixed point of whatever precision makes sense for the application. "Irrational" or "double precision" is irrelevant. 16 bit fixed point should certainly be more than enough, and can be efficiently handled by any microcontroller ever made.
Programs such as Labview of Octave can save you from having to do the actual arithmetic yourself, but they are not a substitute for understanding the math.
-
Do you have values for K1 & K2, and the ranges and precision required for a1, b1, a2 and b2? I would like to play around with them...
Abandoned for now since Labview keeps suggesting irrational numbers for those constants. I tried several times and always got the same result, an this seems to force me to use floating point in double precision, which is really too much for this project, and it's not worth it.
I am going to play with a quadratic approach. It seems more appropriate for my little CPU. Probably the National Instruments curve fitting assistant is made for a PC, where you have floating point in double precision for free as implemented in hardware and with large available resources (lot of ram, high clock frequency, etc).
None of this makes any sense.
Coefficients spat out by a program can be implemented by floating point or fixed point of whatever precision makes sense for the application. "Irrational" or "double precision" is irrelevant. 16 bit fixed point should certainly be more than enough, and can be efficiently handled by any microcontroller ever made.
Programs such as Labview of Octave can save you from having to do the actual arithmetic yourself, but they are not a substitute for understanding the math.
Fully agree.
However, should somebody find themselves in the situation where they do need to calculate fixed-point exponentials please ping me. It isn't that that hard to write a log(x) function using only a lookup table and division, and an exp(x) function using just a table and multiplication.
As pow(a,b) = exp(b*log(a)) this is perfectly achievable in fixed point once the bounds for 'a' and 'b' have been established.
-
However, should somebody find themselves in the situation where they do need to calculate fixed-point exponentials please ping me. It isn't that that hard to write a log(x) function using only a lookup table and division, and an exp(x) function using just a table and multiplication.
The same lookup table in both cases, of exp(2**n), for appropriate +ve and -ve n?
-
As pow(a,b) = exp(b*log(a)) this is perfectly achievable in fixed point once the bounds for 'a' and 'b' have been established.
... also as pow(a, b) = x**(logx(a)*b) for any base x, it's not necessarily best to use e for that base. You can use 10. Or ... if you use 2 as the base then maybe you don't need any table at all (it's just shifts).
-
as pow(a, b) = x**(logx(a)*b) for any base x, it's not necessarily best to use e for that base. You can use 10. Or ... if you use 2 as the base then maybe you don't need any table at all (it's just shifts).
This is a great contribute, and it's probably what I was looking for when I opened this topic! Thanks!
-
For a similar sensor, I was able to manually synthesize an approximation in this terms
distance = (1 / (k1 * Volt + q1))
k1=0.0002391473 (irrational, truncated)
q1=-0.0100251467 (irrational, truncated)
No pow(a,b) required, just a division, a multiplication, and a sum :D
-
Just a proof, but for that similar IR-reflective sensor, this is what I can simplify.
distance = k1 * pow(volt,b1) +q
k1 = 10650.080 (irrational, truncated)
b1 = -0.935 (irrational, truncated)
q = - 10.000
No more three pow, just one.
-
For a similar sensor, I was able to manually synthesize an approximation in this terms
distance = (1 / (k1 * Volt + q1))
k1=0.0002391473 (irrational, truncated)
q1=-0.0100251467 (irrational, truncated)
No pow(a,b) required, just a division, a multiplication, and a sum :D
We're finally starting to get somewhere.
Or try 10000 / (2.391473 * Volt - 100.25)
Or maybe a low order polynomial on the top line also.
-
We're finally starting to get somewhere.
Thanks! Just, don't assume THIS as working model, it's just an example of what I'd like to achieve for the sensor I am using, which is a bit different and it has a different measuring range.
I don't know if it's possible, and that's the reason why I am evaluating other solutions and sensors
I mean, if I apply this approach to my sensor, I get an error of +/-70mm for distances > 1000mm :o :o :o
-
H{ "success", failure" | [ measured=get_value(value), expected ] } :->
(abs(expected - measured) <= tolerance)
H[data, tolerance = 20.00] // 20mm
data{ volt, expected } =
{
{ 2.75, 150.00 }, // volt, mm
{ 2.55, 200.00 },
{ 2.00, 300.00 },
{ 1.55, 400.00 },
{ 1.25, 500.00 },
{ 1.15, 600.00 },
{ 0.90, 700.00 },
{ 0.80, 800.00 },
{ 0.75, 900.00 },
{ 0.65, 1000.00 },
{ 0.60, 1100.00 },
{ 0.55, 1200.00 },
{ 0.50, 1300.00 },
{ 0.44, 1400.00 },
{ 0.45, 1500.00 },
};
get_value(value) :->
k1 = 106500.800
a1 = 1023.000 / 5.000 // 10bit ADC, 5V range
b1 = -0.935
q1 = -100.000
return k1 * pow(a1 * value, b1) + q1
volt=2.750000 failure
expected distance=150.000000
measured distance=185.680132
volt=2.550000 success
expected distance=200.000000
measured distance=206.578036
volt=2.000000 success
expected distance=300.000000
measured distance=284.762786
volt=1.550000 success
expected distance=400.000000
measured distance=388.310394
volt=1.250000 success
expected distance=500.000000
measured distance=497.097488
volt=1.150000 failure
expected distance=600.000000
measured distance=545.510968
volt=0.900000 success
expected distance=700.000000
measured distance=711.781893
volt=0.800000 success
expected distance=800.000000
measured distance=806.289542
volt=0.750000 failure
expected distance=900.000000
measured distance=862.661992
volt=0.650000 success
expected distance=1000.000000
measured distance=1000.479909
volt=0.600000 success
expected distance=1100.000000
measured distance=1085.999996
volt=0.550000 success
expected distance=1200.000000
measured distance=1186.521333
volt=0.500000 success
expected distance=1300.000000
measured distance=1306.433339
volt=0.440000 failure
expected distance=1400.000000
measured distance=1484.994845
volt=0.450000 failure
expected distance=1500.000000
measured distance=1452.038200
A bit better, but ... not so much
-
H[data, tolerance = 60.00] // 60mm
volt=2.750000 success
expected distance=150.000000
measured distance=185.680132
volt=2.550000 success
expected distance=200.000000
measured distance=206.578036
volt=2.000000 success
expected distance=300.000000
measured distance=284.762786
volt=1.550000 success
expected distance=400.000000
measured distance=388.310394
volt=1.250000 success
expected distance=500.000000
measured distance=497.097488
volt=1.150000 success
expected distance=600.000000
measured distance=545.510968
volt=0.900000 success
expected distance=700.000000
measured distance=711.781893
volt=0.800000 success
expected distance=800.000000
measured distance=806.289542
volt=0.750000 success
expected distance=900.000000
measured distance=862.661992
volt=0.650000 success
expected distance=1000.000000
measured distance=1000.479909
volt=0.600000 success
expected distance=1100.000000
measured distance=1085.999996
volt=0.550000 success
expected distance=1200.000000
measured distance=1186.521333
volt=0.500000 success
expected distance=1300.000000
measured distance=1306.433339
volt=0.440000 failure
expected distance=1400.000000
measured distance=1484.994845
volt=0.450000 success
expected distance=1500.000000
measured distance=1452.038200
"expected" means the value in the table.
"measured" means the compute value that should approximate it.
-
.
-
Yup, this topic is a solid amount of embedded things ... hacking of things, inner digital math on "grinder old CPU", and crazy stuff :D
find a way to implement pow() in a microcontroller not using the supplied one
find a way to fit your data to a given expression (adjusting some parameters)
find an expression that fits your data
{ Yes, Yes, Yes }
Plus I am open to alternative solutions! For example, yesterday it was suggested to use a different approach for the sensor
Have you considered an optical encoder -- either linear or rotary?
and now I'm also thinking of opening the motor box and hacking it to get a quadrature sensor. I am very reluctant to do this because it is a very dirty trick, the motor box is sealed and I need to "Dremel cut", and, since a quadrature encoder is a relative-position encoder it would require the frame to "Go back home" to the "zero position" (600mm from the ground) in order to set the zero, which is also something I'd like to avoid.
Go Back home and set the zero -> offset = 600mm
now start counting, position = offset + encoder'delta * counted-steps
Dunno, I will try it the next weekend :o
-
Here's something really basic, that calculate log() and exp() base 2 over limited range using a really silly mix of double and floating point numbers, but should give ideas of how you could implement them when you know what range you will be working over.
#include <stdio.h>
double exp_table [16] = {
1.000169239705302138, 1.000338508052682318, 1.000677130693066408, 1.001354719892108225,
1.002711275050202522, 1.005429901112802726, 1.010889286051700475, 1.021897148654116627,
1.044273782427413755, 1.090507732665257690, 1.189207115002721027, 1.414213562373095145,
2.000000000000000000, 4.000000000000000000, 16.00000000000000000, 256.0000000000000000,
};
double my_pow2(double a) {
int fixed = a * 4096;
double total = 1.0;
for(int i = 0; i < 16; i++) {
if(fixed & (1<<i))
total *= exp_table[i];
}
return total; // 2^(a/4096) for a limited range
}
double log_table [16] = {
1.00000528830729185, 1.00002115339696473, 1.00004230724139576, 1.00008461627269440,
1.00016923970530214, 1.00033850805268232, 1.00067713069306641, 1.00135471989210822,
1.00271127505020252, 1.00542990111280273, 1.01088928605170048, 1.02189714865411663,
1.04427378242741375, 1.09050773266525769, 1.18920711500272103, 1.41421356237309515
};
int my_log2_fixed(double a) {
int fixed = 0;
while(a >= 2.0) {
fixed++;
a /= 2.0;
}
while(a < 1.0) {
fixed --;
a *= 2.0;
}
fixed <<=16;
for(int i = 15; i >= 0; i--) {
if(a > log_table[i]) {
a /= log_table[i];
fixed |= 1<<i;
}
}
return fixed; // log2(a) * 65536
}
-
But why bother curve-fitting at all - just interpolate between points.
#include <stdio.h>
struct {
double x, y, k;
} data[] =
{
{ 2.75, 150.00 }, // volt, mm
{ 2.55, 200.00 },
{ 2.00, 300.00 },
{ 1.55, 400.00 },
{ 1.25, 500.00 },
{ 1.15, 600.00 },
{ 0.90, 700.00 },
{ 0.80, 800.00 },
{ 0.75, 900.00 },
{ 0.65, 1000.00 },
{ 0.60, 1100.00 },
{ 0.55, 1200.00 },
{ 0.50, 1300.00 },
{ 0.44, 1400.00 }
};
double get_value(double i_value) {
int i;
// Find where we are interpolating
for(i = 0; i < sizeof(data)/sizeof(data[0])-2; i++) {
if(i_value > data[i+1].x) {
break;
}
}
return data[i].y + (i_value-data[i].x) * data[i].k;
}
int main(int argc, char *argv[]) {
// Calc some constants
for(int i = 0; i < sizeof(data)/sizeof(data[0])-1; i++) {
data[i].k = (data[i+1].y-data[i].y) / (data[i+1].x-data[i].x);
}
// Check against the table values
for(int i = 0; i < sizeof(data)/sizeof(data[0]); i++) {
printf("%5.2f, %7.2f, %7.2f\n", data[i].x, data[i].y, get_value(data[i].x));
}
// Check between two numbers
printf("\n%5.2f, %7.2f, %7.2f\n", 0.85, 750.00, get_value(0.85));
}
-
I accidentally searched:
https://html.duckduckgo.com/html?q=fractional+look+up+table (https://html.duckduckgo.com/html?q=fractional+look+up+table)
and found:
https://www.sciencedirect.com/topics/computer-science/fractional-value (https://www.sciencedirect.com/topics/computer-science/fractional-value)
I like the idea of first scaling a number between 0.5 and 1 and then using a LUT or simplified approximation.
But scaling it back would need an exponential function so I'm confuses as to whether this is applicable in this situation.
-
Why not simply approximate the transfer function with a polynomial, as already suggested by rstofer?
A 5th order polynomial fits nicely (see plot). I would not use a higher order as the provided data are ovbiously noisy.
The evaluation requires 5 multiplications and 6 additions then. The coefficients can be either be stored in a table, or immediate values in the code.
>> data
data =
2.75000 15.00000
2.55000 20.00000
2.00000 30.00000
1.55000 40.00000
1.25000 50.00000
1.15000 60.00000
0.90000 70.00000
0.80000 80.00000
0.75000 90.00000
0.65000 100.00000
0.60000 110.00000
0.55000 120.00000
0.50000 130.00000
0.44000 140.00000
0.45000 150.00000
>> P = polyfit(data(:,1), data(:,2), 5)
P =
-12.951 118.581 -429.858 782.454 -750.827 357.198
>> x = 0.44:0.01:2.75;
>> plot(data(:,1), data(:,2), ".", x, polyval(P,x)); grid on
-
Why not simply approximate the transfer function with a polynomial, as already suggested by rstofer?
Good idea! I am still studying Octave, and I did a couple of attempts, but then I have been busy with a complete different job task.
provided data are obviously noisy.
The nature of the sensor is IR reflective, that's very noisy by itself, that's also the reason why I want to filter and to amplify the signal. Not yet done, but I will this weekend.
I am also waiting for a couple of better ADC at 12 and 16 bit of resolution. I need to write an SPI function to manage them, but hey? look at the current firmware only takes 512 byte, and it manages a mean-window filtering function, the ADC conversion, the serial talking to the host (with a minimal protocol to setup-up working variables), and a 100uS-fine timer used to control the motor.
All written in assembly. I love it!
I am going to reuse the 100uS-fine timer for a new ultrasonic ranging sensor used with quad rotors. It's modulator is smarter than traditional ultrasonic ranging sensors, and its outputs is a digital pulse of the width directly proportional to the measured distance.
I ordered it two days ago, it will probably arrive tomorrow.
So I may use both: IR reflective + ultrasonic, or just one. We will see :D
-
I often use Gnuplot for fitting data to a function I define myself.
For example, if you have data.txt with input value in the first column, output value in the second column, and if there are any comments they start with a #, you can start Gnuplot, and tell it e.g.
f(x) = C0 + x*(C1 + x*(C2 + x*(C3 + x*(C4 + x*C5))));
fit f(x) "data.txt" using 1:2 via C0,C1,C2,C3,C4,C5
and you can also draw both the sample data and the curve fit using e.g.
plot "data.txt" using 1:2 notitle with points lc "black", f(x) notitle with lines lc "red"
You can print the coefficients using e.g.
print(C0, C1, C2, C3, C4, C5)
but Gnuplot will automatically save a text file called fit.log which contains the fitted variables (C0 through C5, above) under the "Final set of parameters" heading; and all kinds of convergence, error, and correlation matrix information that you may or may not find useful.
If you want the function in a variable name to be used in the plot, you can use sprintf (which has the same syntax the C printf() function has):
fstr = sprintf("f(x) = %f + x*(%f + x*(%f + x*(%f + x*(%f + x*%f))))", C0, C1, C2, C3, C4, C5)
so that
set title fstr ; plot "data.txt" using 1:2 notitle with points lc "black", f(x) notitle with lines lc "red"
will add the function definition as the title of the plot.
I've found this to be very useful when I am still exploring what kind of a function fits my needs most. (In Gnuplot, AB is written A**B, and you can use conditional expressions – (condition ? value-if-true : value-if-false) – in functions. ^|& are binary operators between integers (exclusive-or, or, and, respectively), so 2^5==7 and not 32. == tests for equalty, and != for inequality.)
As an example, when dealing with distances r, I often try to avoid r and 1/r, preferring r2 and 1/r2 instead, because if r is calculated from a vector, odd powers need a "slow" square root, whereas even powers are "cheap" (just need multiplications; and only one division for all even powers of the reciprocal of the distance).
When implementing the function in fixed point format, int(expression) is very useful, because it truncates the expression to an integer. However, the fit variables (C0 through C5, above) must still be treated as reals, as otherwise Gnuplot cannot converge to a solution. (This means, if you try to use int(C0) in f(x), you can expect Gnuplot to not be able to converge to a solution, because small changes in C0 do not have any effect on the result.)
-
For most of my "manual" (meaning not automatically generated - for this Gnuplot is nice, there are others) graphing needs, I use Scidavis. Many features including curve fitting.
-
automatically generated - for this Gnuplot is nice
Indeed, the way I usually store the derivation of the coefficients is to write the fitting (and plotting) into a script. But mostly, it just fits my workflow.
I might be strange, but I don't use Jupyter notebooks (for Sagemath), wxMaxima documents, Octave workspaces, Maple notebooks, etc.; I prefer to save the operations into a text file. It is faster for me to less or find+grep details to find the particular project, than to open a dozen or so of them in their respective programs to see the content. Then again, I have hundreds of those snippets...
For proper documentation, there's LibreOffice and LaTeX and other, better tools.
-
I have been very busy with other stuff but I haven't forgotten this project.
Today I managed to find two hours to modify my board and to implement an new uploader software able to auto-reset the target. The RS232 signal CTS is used for this.
uploader-v1 # ./obj/binvect-uart-upld hAllo.bin
loading binvect ... done, 512 bytes loaded
opening /dev/ttyS0 at 8N1-1200bps ... success
resetting target and waiting for ack ... done
uploading: [########################################] 100%, success
It works, and it allows me to develop things remotely and to upload them to the target so I can work on it, collects data, experiment stuff, even when I am outdoors :D
This uploader also works with the chip in expanded mode. I have already prepared a new board and a special "first stage bootloader" and it's ready to work.
If I don't have to go overseas for work, we will see some progress soon.