C Math wiz required please...

[liquidsquid]

Hi to any math wizards out there!

I am trying to find an efficient way to calculate LOG10(), POW() (only POW(10,?)) and SQRT(). Or, use tables... my program is far too slow on an embedded controller that has no floating point core.

I need to step a log value by a linear one, or in other words, step RF levels in dBm by volts.


The formula to convert volts to dBm is:

dBm = 10 * log(((V^2)/50ohm)/1mW)

Or simplified: dBm = 10 * log10(20 * V * V)


and dBm to volts:

volts = 0.223607 * sqrt(pow(10, dBm / 10))

Where 0.223.. = Volts at 1mW

Now dBm I have multiplied by 10,000 to use integer math rather than double-precision floats which speeds up text conversions considerably for display. For these routines I have to convert back to double for the built-in C routines, which is where my time gets chewed up, almost 100mS to step once!

I have found some good, fast sqrt() integer routines at http://www.azillionmonkeys.com which I have used in other programs before like DFTs and FFTs, but was pointless here until I figured out how to get all results into integer form with decent precision.

Currently this is the long-hand version of code (just a portion of a much larger multiple units conversion routine):

Code: Select all

 result = dBm10ktoVolts(level->dBm10k);      //Convert to volts
 result += levelStepV;                     //Step by volts
 level->dBm10k = VoltstodBm10k(result);  //Convert back to dBm10k


Code: Select all

/************ Level conversion routines, common base of dBm*1000 *****************/
/******* dB are in the format of -000.0000dB, so a fixed dP at place 4 **/
/******* uV are in the format of 0uV, so no dP or 000.000000V**/
/******* uW are in the format of 000,000,000uW, so no dP, or 000.000000W ****/

double dBm10ktoVolts(INT32 dBm10k){
   return(0.22360679774997896964 * sqrt(pow(10.0, dBm10k / 100000.0)));
}

INT32 VoltstodBm10k(double Volts){
   return(100000.0 * log10(20.0 * Volts * Volts)); //Convert to dBm * 1000 (10 * 1000 for the 10000)
}


Anyone have any linkage for integer-based log10 or pow routines? I could go from there. Otherwise, any clever ideas I can get started on?

Log10 can be easy in integer math since the value repeats between powers of 10, and you just need a table lookup and some interpolation to get close... but Exp(10,?) seems to have my poor brain befuddled.

Thanks for some ideas...

See what old PC programmers had to deal with in SX/DX days? Poor guys...

-LS

[liquidsquid]

Oh, wanted to mention that I was an idiot with math, as I should have made everything dBm32K rather than 10K so I could use shifts rather than /10.

Fixing as I type...

[IntelMole]

Just so I know what you're going on about:

You're after a routine that will calculate the logs of a given number.

My first instinct says power series, then divide by ln(10), or 2.3026 (5sf) to get into the correct base. The only thing that might screw that is the radius of convergance for the power series of a logarithm... Anyone know it?

If it did converge, it's only accurate to how many terms you use, which depends on how much time you have to calculate it...

I did know at one point two weeks ago how to do power series (when I had it for a revision question in preparation for an exam) so if you really wanna do it that way, then I'll dig out.

As for the power thing, I'm not to sure what the fastest method is, maybe I'll have an inspiration somewhere ,
-Mole

[Flowboy]

You can simplify your expression some more by using the identity that log(ab) = log(a) + log(b)

So you have:

dBM = 10*log(20*V*V)

but this will simplify to:

dBM = 10*(log(20) + 2*log(v))

You've now removed any need for power or square root.

My instinct for an embedded system would be to change everything to log base 2.

Take a look at this site for a good discussion of log and ways of approximating it.

http://www.dattalo.com/technical/theory/logs.html

[UberGerbil]

Hi to any math wizards out there!
See what old PC programmers had to deal with in SX/DX days? Poor guys...
-LS

Well, remember that prior to the 486, nobody had an x87 except for the hardcore who had actually gone out and bought the separate coprocessor. So mostly we used math packages that either emulated the x87 (and used it the machine actually happened to have it available) or used the alternate math library that did it all in software (using different algorithms that were faster than the emulated x87 but slower than the actual hardware). The MS C compiler used to offer those choices under a switch; I haven't looked recently so I don't know if it still does.

You may be able to find the equivalent libraries... something like this.

So even back in the good 'ol days we mostly didn't have to translate floating point code into pure integer ops. Though I do remember writing a fractal program entirely using ints (with the inner loop in assembly).

Flowboy has some good suggestions. Log2 will be better if you can use it. Try to avoid division by all means -- precalculate recipricals and use that instead. If you're not going to be dealing with negative values you can gain an extra bit of precision using UINTs.

[Asin]

Once you have the pow() method/function, couldn't you just use Taylor series and Taylor approximations using derivates for the rest?

[Flowboy]

pow(10, x) is doable with a lookup table too.

Use the following method:

Example:

10^42 = 10^40*10^2

10^357 = 10^300*10^50*10^7

So what you do is you do a divide by 100, and get the quotient and the remainder. Use the quotient as an index into a lookup table of precalculated values (10^100, 10^200, 10^300....). So for our example we have a quotient of 3 which gives us 10^300 as the precalculated value.

Now take the remainder, in our example, 57. Do a divide by 10 and get the quotient and remainder. Again, take the quotient as an index into a lookup table, which gives us 10^50.

Take the remainder, which should be 7, and use this as an index into a lookup table. We get 10^7.

Now multiply those guys together. We get 10^300*10^50*10^7 = 10^357

Now here's the clever part. How about we do it in base 2?

Now let's see what happens when we use powers of 2 as the number to divide and get the remainder by instead of 10.

I'm going to use a smaller number to simplify working. We want to find 10^9.

8 goes into 9 once with a remainder 1. So we have 10^8 to remember.
4 goes into 1 zero times with a remainder 1. So we have no result this pass.
2 goes into 1 zero times with a remainder 1. So we have no result this pass.
1 goes into 1 one time with no remainder. So we have 10^1 this time round.

So we have 10^8*10^1 as the result, which = 10^9

Notice something?

We can just test the bits as we go.

Code: Select all

int nExponent; // your input
int nBitmask = 1;
int nResult = 1;
int nPowerOf10 = 10;
for (int i = 0; i < 16; i++)
{
if (nExponent & nBitmask)
{
nResult *= nPowerOf10;
}

nBitmask <<= 1;
nPowerOf10 *= 10;
}


[the_silver_bullet24]

C Math? Well I guess the C doesn't stand for consumer

[IntelMole]

If you're after a power series of logs:
http://mathworld.wolfram.com/NaturalLogarithm.html

log2(1+x) = (x - .5x + .33x^2 - .25x^4 + .2x^5) / ln 2

or...

log2(1+x) = (x - .5x + .33x^2 - .25x^4 + .2x^5) / (0.693...)

Enjoy coding that ,
-Mole

[Flowboy]

Liquidsquid,

Can you give me some idea of the dynamic range of the voltage and dbM measurements? Maximum and minimum values and required precision. Presumably these are limited by the ADC / DAC components you're using.

We might be able to do some lookup table fun.

[just brew it!]

How many steps are you talking about here, and how constrained are you by the amount of memory available?

At first glance, it sounds like something that should be done using a pre-computed lookup table. If you need a lot of steps, and memory is really tight, you might be able to get away with storing only every 2nd or 3rd point and doing a linear interpolation (depending on how accurate you need to be).

[liquidsquid]

Ah, I haven't been back to look at this thread in a bit, a lot of interesting ideas. I was working on this the past Friday (my first opportunity) and have decided on the table lookup method. The range of values is from -140 dBm to 0dBm

The formula of dBm to volts of:

0.22360679774997896964 * sqrt(pow(10.0, dBm10k / 100000.0))

In actual non-scaled format is:

volts = 0.22360679774997896964 * sqrt(10^(dBm/10))

The 0.223... is a 1mW reference point.

Anyhow, foretting my math completely, I realized that

sqrt(10^x) == 10^(x/2) since sqrt(x) == x^(1/2)

So I can skip the sqare root.

So I take my dBm value, divide it by 20 rather than 10. Instant sqare root.

volts = 0.22360679774997896964 * 10^(dBm/20)

Now I take the value in a long int (32 bits in my case) and set it in 16.16 fixed point format (meaning I shift the value 16 places or multiply by 2^16). Then the long int is unionized with two unsigned words so it is easy to extract the integer portion and the fractional just by refering to each portion of the long int. Each portion can be used to refer to a lookup table. The int portion being the scaler (power of 10), and the fractional, the log table, which repeats for each increment of the scaler. The fractional lookup table has 256 values to choose from using the upper byte of the fractional, and then it is further calculated by interpolating between the two points it is closest to, with the second byte of the fractional to weight each value (pretty close, but not perfect, but I don't need perfection).

Then I simply divide the value returned by the fractional portion by the value retrieved by the int portion (pre-multiplied by 1/0.223...), and I have my result. I will post the source later this week. I went from 45mS per calculation to <100uS, though I haven't accuratly benchmarked it yet.

The accuracy however is limited by the smallest value a 16.16 format number can hold, or 1/(2^16) or 15uV, so I have to play with the scaler a bit more, or return the value as a float to increase accuracy since -140dBm is much smaller than this. I will take a hit in the float conversion, but nothing close to the hit as using the POW() function.

-LS

[liquidsquid]

I must have been tired when I typed the original post, the actual formula is:

Volts = V1mW * 10^(dBm/20) after simplification, and V1mW = 0.223606...

Much easier that what I posted (I changed the post so the original f-up isn't there).

This is the result for the dBm to float. Just a bit faster . Just started on the other direction, which should be similar.

Code: Select all

//For the 10^x in 256 steps, in reverse order since we work in negatives, it saves a subtraction
//step.
const float powTableFloat[257] = {10.0,9.9104585625,9.8217188919,9.7337738090,9.6466161991,9.5602390110,
   9.4746352566,9.3897980105,9.3057204093,9.2223956510,9.1398169947,9.0579777594,8.9768713245,
   8.8964911282,8.8168306678,8.7378834985,8.6596432336,8.5821035433,8.5052581544,8.4291008503,
   8.3536254696,8.2788259063,8.2046961090,8.1312300806,8.0584218776,7.9862656097,7.9147554394,
   7.8438855815,7.7736503024,7.7040439201,7.6350608034,7.5666953714,7.4989420933,7.4317954878,
   7.3652501227,7.2993006144,7.2339416274,7.1691678741,7.1049741144,7.0413551549,6.9783058486,
   6.9158210949,6.8538958387,6.7925250701,6.7317038241,6.6714271804,6.6116902624,6.5524882374,
   6.4938163158,6.4356697510,6.3780438389,6.3209339175,6.2643353666,6.2082436072,6.1526541015,
   6.0975623522,6.0429639024,5.9888543349,5.9352292723,5.8820843762,5.8294153471,5.7772179241,
   5.7254878844,5.6742210428,5.6234132519,5.5730604013,5.5231584173,5.4737032629,5.4246909370,
   5.3761174746,5.3279789459,5.2802714565,5.2329911468,5.1861341918,5.1396968008,5.0936752168,
   5.0480657167,5.0028646106,4.9580682417,4.9136729859,4.8696752517,4.8260714794,4.7828581417,
   4.7400317423,4.6975888167,4.6555259312,4.6138396827,4.5725266990,4.5315836376,4.4910071863,
   4.4507940624,4.4109410125,4.3714448126,4.3323022674,4.2935102101,4.2550655025,4.2169650343,
   4.1792057232,4.1417845144,4.1046983804,4.0679443211,4.0315193629,3.9954205589,3.9596449889,
   3.9241897585,3.8890519993,3.8542288686,3.8197175493,3.7855152493,3.7516192015,3.7180266639,
   3.6847349187,3.6517412725,3.6190430563,3.5866376245,3.5545223556,3.5226946515,3.4911519372,
   3.4598916609,3.4289112936,3.3982083289,3.3677802831,3.3376246943,3.3077391230,3.2781211514,
   3.2487683834,3.2196784443,3.1908489806,3.1622776602,3.1339621714,3.1059002236,3.0780895465,
   3.0505278903,3.0232130250,2.9961427410,2.9693148482,2.9427271762,2.9163775741,2.8902639100,
   2.8643840715,2.8387359648,2.8133175149,2.7881266654,2.7631613785,2.7384196343,2.7138994312,
   2.6895987856,2.6655157314,2.6416483204,2.6179946216,2.5945527214,2.5713207234,2.5482967480,
   2.5254789326,2.5028654312,2.4804544143,2.4582440689,2.4362325982,2.4144182213,2.3927991734,
   2.3713737057,2.3501400846,2.3290965925,2.3082415268,2.2875732003,2.2670899410,2.2467900918,
   2.2266720104,2.2067340691,2.1869746550,2.1673921696,2.1479850285,2.1287516618,2.1096905134,
   2.0908000413,2.0720787172,2.0535250265,2.0351374682,2.0169145547,1.9988548119,1.9809567786,
   1.9632190068,1.9456400616,1.9282185208,1.9109529750,1.8938420273,1.8768842936,1.8600784018,
   1.8434229924,1.8269167179,1.8105582430,1.7943462442,1.7782794100,1.7623564406,1.7465760477,
   1.7309369547,1.7154378963,1.7000776188,1.6848548794,1.6697684466,1.6548170999,1.6399996297,
   1.6253148373,1.6107615346,1.5963385443,1.5820446995,1.5678788438,1.5538398313,1.5399265261,
   1.5261378026,1.5124725453,1.4989296487,1.4855080172,1.4722065648,1.4590242156,1.4459599031,
   1.4330125702,1.4201811697,1.4074646633,1.3948620224,1.3823722274,1.3699942677,1.3577271421,
   1.3455698581,1.3335214322,1.3215808896,1.3097472643,1.2980195990,1.2863969449,1.2748783618,
   1.2634629177,1.2521496891,1.2409377608,1.2298262257,1.2188141848,1.2079007474,1.1970850305,
   1.1863661591,1.1757432659,1.1652154917,1.1547819847,1.1444419008,1.1341944035,1.1240386638,
   1.1139738600,1.1039991779,1.0941138106,1.0843169582,1.0746078283,1.0649856354,1.0554496009,
   1.0459989534,1.0366329284,1.0273507682,1.0181517217,1.0090350448,1.0000000000};

//Enough values for -170dBm, shifted + 1 since normal values are negative, but
//We look up in the positive realm. (4.472135955 is tossed representing 0)
const float powScalersFloat[16] = {44.72135955,447.2135955,4472.1359549996,
   44721.3595499958,447213.595499958,4472135.954999579,44721359.5499958,447213595.499958,
   4472135954.999579,4.4721359550e10,4.4721359550e11,4.4721359550e12,4.4721359550e13,
   4.4721359550e14,4.4721359550e15};

/***********************************************************************************
** Converts from dBm65k (16.16 fixed point int) to volts as a float.
** Simplifies the formula of: Volts := .22360679774997896964 * sqrt(10^(dBm/10))
** to two divides (one int, one float) two lookups.
** only handles the range from 0 to -160dBm, not positives!
***********************************************************************************/   
float fastPow10dBm(union dBm65kUnion dBm65k) {
   
   float unscaled;
   float scaler;
   
   dBm65k.dBm /= -20;                  //Formula is sqrt(10^(x/10)) so divide by twenty, and convert to positive.
                                 //sqrt(10^(x/10)) == 10^(x/20)
      
   unscaled = powTableFloat[dBm65k.bytes.largerfrac];      //Ignore the lower fractional value during lookup
                                 //We can increase accuracy (lowering poerformance) by
                                 //Interpolating between the lookup value and the previous
                                 //value in the table. No need found yet...
   scaler = powScalersFloat[dBm65k.words.integer];   
      
   return unscaled / scaler;
}


[liquidsquid]

Oh, and the volts to dBm also works out to be:

= 20 * log(Vin / V1mW), where the V1mw = 0.223....

So

20 * (log(Vin) - log(V1mW))

The only portion I have to deal with then is the log(Vin), which will be fun (not!) since it has to handle quite a dynamic range. Perhaps I can yank apart the float using a union to tap into the portions of the floats. So far so good...

-LS

[just brew it!]

For log(Vin) maybe you could use a table containing input/output value pairs, and do a binary search.

[liquidsquid]

Got it, this is the other direction. I have to increase precision however, so I will take a hit during interpolation.

Code: Select all

//This is the values for the numeric range of log(0.5) to log(1.0), which for every value
//This gives a range of 0-1 (not all values are scaled by 20 * 2^16 / (log(10,2))
//All set with log(V1mw) already added in.
const long int log2Table_p5_to_1[257] = {458077,460296,462507,464709,466903,469088,471265,473433,475593,
   477746,479890,482026,484154,486274,488386,490490,492587,494676,496757,498831,500897,502956,505007,
   507051,509088,511117,513139,515154,517162,519163,521157,523144,525124,527097,529063,531023,532975,
   534921,536861,538794,540720,542640,544553,546460,548361,550255,552143,554025,555901,557770,559633,
   561491,563342,565187,567026,568860,570687,572509,574324,576134,577939,579737,581530,583317,585099,
   586875,588646,590411,592170,593924,595673,597417,599155,600888,602615,604338,606055,607767,609474,
   611175,612872,614564,616250,617932,619609,621281,622947,624609,626267,627919,629567,631209,632847,
   634481,636110,637734,639353,640968,642578,644184,645785,647382,648974,650562,652146,653725,655299,
   656870,658436,659997,661555,663108,664657,666201,667742,669278,670811,672339,673863,675383,676899,
   678411,679919,681423,682923,684419,685911,687399,688883,690364,691840,693313,694782,696247,697709,
   699167,700621,702071,703518,704960,706400,707835,709267,710696,712121,713542,714960,716374,717785,
   719192,720596,721996,723393,724787,726177,727563,728947,730327,731703,733077,734447,735813,737177,
   738537,739894,741248,742598,743946,745290,746631,747968,749303,750635,751963,753288,754611,755930,
   757246,758559,759869,761176,762481,763782,765080,766375,767667,768957,770243,771527,772807,774085,
   775360,776632,777901,779167,780431,781692,782950,784205,785457,786707,787954,789198,790440,791679,
   792915,794148,795379,796607,797833,799055,800276,801493,802708,803921,805131,806338,807543,808745,
   809945,811142,812336,813528,814718,815905,817090,818272,819452,820629,821804,822977,824147,825314,
   826480,827643,828803,829961,831117,832270,833422,834570,835717,836861,838003,839143,840280,841415,
   842548,843678,844807,845933,847057,848178,849298,850415,851530,852643};
   
/***********************************************************************************
** Convert volts to dBm65k format (16.16)
** formula is dBm := 20 * log(Vin / V1mw)
** Use the rules log(ab) = log(a) + log(b)
** and log(a/b) = log(a) - log(b)
** so dBm = 20 * (log(Vin) - log(V1mW))
** volts are always positive.
***********************************************************************************/
union dBm65kUnion fastVtodBm65k(double volts) {
   union dBm65kUnion result;
   unsigned long int normalized;
   int exponent2;
   long int kickass;
      
   normalized = (unsigned long int)(ldexp(frexp(volts, &exponent2), 9)) - 256;   //Normalized = 0-256, or 1/2 to 1. Int is power of 2.

   kickass = exponent2 * 394566;         //Scaler for 2^16 * 20 * log(Vin, 2) / log(10, 2)
                                 //394566 = 2^16 * 20 / log(10, 2)                              
   if (normalized > 0) {
      kickass += log2Table_p5_to_1[normalized];
   } else {
      kicakass = 0x7FFFFFFF;            //Huge, doesn't really matter
   }
   result.dBm = kickass;
   return result;
}


-Ls

[liquidsquid]

Oh, btw I re-wrote these as I had some rounding errors to deal with and precision problems. They are more general now and reusable, so I will post when I get back to work.

-LS