0

I've implemented linear interpolation on my PIC18F device where I'm interpolating ADC values with intervals stored in EEPROM. I tested the code and it doesn't calculate properly. So I put it into my computer and the calculations works well so I don't really know where the problem might come from, as the two functions I use right now for testing purposes are same on both platforms.

These are the functions same for both platforms:

uint16_t interpolate(uint16_t a, uint16_t b, uint16_t c, uint16_t d, uint16_t raw) {
    // b+(x-a)*(d-b)/(c-a)
    return b + (raw - a) * (d - b) / (c - a);
}

uint16_t mapRawToCal(uint16_t raw) {
    uint16_t cal = 0;

    if (raw < 150) {
        cal = 0;
    }
    else if (raw >= 150 && raw <= 457) {
        cal = interpolate(150, 0, 457, 800, raw);
    }
    else if (raw > 457 && raw <= 1235) {
        cal = interpolate(457, 800, 1235, 1000, raw);
    }
    else if (raw > 1235 && raw <= 1988) {
        cal = interpolate(1235, 1000, 1988, 1988, raw);
    }
    else if (raw > 1988 && raw <= 2700) {
        cal = interpolate(1988, 1988, 2700, 3500, raw);
    }
    else if (raw > 2700 && raw <= 3560) {
        cal = interpolate(2700, 3500, 3560, 3300, raw);
    }
    else if (raw > 3560 && raw <= 3987) {
        cal = interpolate(3560, 3300, 3987, 4096, raw);
    }
    else if (raw > 3987) {
        cal = 4096;
    }

    return cal;
}

This is the PIC18F device main function:

int main(void) {
    // Setup
 hardware_init();
   
 ADC_Init();
    
    // Loop
    while(1) {
        uint16_t ch1Data;
        ch1Data = ADC_Read(1);
        
        uint16_t *buf = usb_get_in_buffer(1);
        buf[0] = ch1Data; // Rx Axis
        buf[1] = mapRawToCal(ch1Data); // Ry Axis
           
        usb_send_in_buffer(1, 4);
           
        __delay_ms(200); //Delay

        #ifndef USB_USE_INTERRUPTS
        usb_service();
        #endif
    }
    
 return 0;
}

And this is the main function on Windows:

int main() {
    for (uint16_t raw = 0; raw <= 4096; raw++) {
        printf("input=%d - output=%d\n", raw, mapRawToCal(raw));
    }

    return 0;
}

Here is chunk of the C Windows program result in console (that is accurate):

input=129 - output=0
input=130 - output=0
input=131 - output=0
input=132 - output=0
input=133 - output=0
input=134 - output=0
input=135 - output=0
input=136 - output=0
input=137 - output=0
input=138 - output=0
input=139 - output=0
input=140 - output=0
input=141 - output=0
input=142 - output=0
input=143 - output=0
input=144 - output=0
input=145 - output=0
input=146 - output=0
input=147 - output=0
input=148 - output=0
input=149 - output=0
input=150 - output=0
input=151 - output=2
input=152 - output=5
input=153 - output=7
input=154 - output=10
input=155 - output=13
input=156 - output=15
input=157 - output=18
input=158 - output=20
input=159 - output=23
input=160 - output=26
input=161 - output=28
input=162 - output=31
input=163 - output=33
input=164 - output=36
input=165 - output=39
input=166 - output=41
input=167 - output=44
input=168 - output=46
input=169 - output=49
input=170 - output=52
input=171 - output=54
input=172 - output=57
input=173 - output=59
input=174 - output=62
input=175 - output=65
input=176 - output=67
input=177 - output=70
input=178 - output=72
input=179 - output=75
input=180 - output=78
input=181 - output=80
input=182 - output=83
input=183 - output=85
input=184 - output=88
input=185 - output=91
input=186 - output=93
input=187 - output=96
input=188 - output=99
input=189 - output=101
input=190 - output=104
input=191 - output=106
input=192 - output=109
input=193 - output=112
input=194 - output=114
input=195 - output=117
input=196 - output=119
input=197 - output=122
input=198 - output=125
input=199 - output=127
input=200 - output=130
input=201 - output=132
input=202 - output=135
input=203 - output=138
input=204 - output=140
input=205 - output=143
input=206 - output=145
input=207 - output=148
input=208 - output=151
input=209 - output=153
input=210 - output=156
input=211 - output=158
input=212 - output=161
input=213 - output=164
input=214 - output=166
input=215 - output=169
input=216 - output=171
input=217 - output=174
input=218 - output=177
input=219 - output=179
input=220 - output=182
input=221 - output=185
input=222 - output=187
input=223 - output=190
input=224 - output=192
input=225 - output=195
input=226 - output=198
input=227 - output=200
input=228 - output=203
input=229 - output=205
input=230 - output=208
input=231 - output=211
input=232 - output=213
input=233 - output=216
input=234 - output=218
input=235 - output=221
input=236 - output=224
input=237 - output=226
input=238 - output=229
input=239 - output=231
input=240 - output=234
input=241 - output=237
input=242 - output=239
input=243 - output=242
input=244 - output=244
input=245 - output=247
input=246 - output=250
input=247 - output=252
input=248 - output=255
input=249 - output=257
input=250 - output=260
input=251 - output=263
input=252 - output=265
input=253 - output=268
input=254 - output=271
input=255 - output=273
input=256 - output=276
input=257 - output=278
input=258 - output=281
input=259 - output=284
input=260 - output=286
input=261 - output=289
input=262 - output=291
input=263 - output=294
input=264 - output=297
input=265 - output=299
input=266 - output=302
input=267 - output=304
input=268 - output=307
input=269 - output=310
input=270 - output=312
input=271 - output=315
input=272 - output=317
input=273 - output=320
input=274 - output=323
input=275 - output=325
input=276 - output=328
input=277 - output=330
input=278 - output=333
input=279 - output=336
input=280 - output=338
input=281 - output=341
input=282 - output=343
input=283 - output=346
input=284 - output=349
input=285 - output=351
input=286 - output=354
input=287 - output=357
input=288 - output=359
input=289 - output=362
input=290 - output=364
input=291 - output=367
input=292 - output=370
input=293 - output=372
input=294 - output=375
input=295 - output=377
input=296 - output=380
input=297 - output=383
input=298 - output=385
input=299 - output=388
input=300 - output=390
input=301 - output=393
input=302 - output=396
input=303 - output=398
input=304 - output=401
input=305 - output=403
input=306 - output=406
input=307 - output=409
input=308 - output=411
input=309 - output=414
input=310 - output=416
input=311 - output=419
input=312 - output=422
input=313 - output=424
input=314 - output=427
input=315 - output=429
input=316 - output=432
input=317 - output=435
input=318 - output=437
input=319 - output=440
input=320 - output=442
input=321 - output=445
input=322 - output=448
input=323 - output=450
input=324 - output=453
input=325 - output=456
input=326 - output=458
input=327 - output=461
input=328 - output=463
input=329 - output=466
input=330 - output=469
input=331 - output=471
input=332 - output=474
input=333 - output=476
input=334 - output=479
input=335 - output=482
input=336 - output=484
input=337 - output=487
input=338 - output=489
input=339 - output=492
input=340 - output=495
input=341 - output=497
input=342 - output=500
input=343 - output=502
input=344 - output=505
input=345 - output=508
input=346 - output=510
input=347 - output=513
input=348 - output=515
input=349 - output=518
input=350 - output=521
input=351 - output=523
input=352 - output=526
input=353 - output=528
input=354 - output=531
input=355 - output=534
input=356 - output=536
input=357 - output=539
input=358 - output=542
input=359 - output=544
input=360 - output=547
input=361 - output=549
input=362 - output=552
input=363 - output=555
input=364 - output=557
input=365 - output=560
input=366 - output=562
input=367 - output=565
input=368 - output=568
input=369 - output=570
input=370 - output=573
input=371 - output=575
input=372 - output=578
input=373 - output=581
input=374 - output=583
input=375 - output=586
input=376 - output=588
input=377 - output=591
input=378 - output=594
input=379 - output=596
input=380 - output=599
input=381 - output=601
input=382 - output=604
input=383 - output=607
input=384 - output=609
input=385 - output=612
input=386 - output=614
input=387 - output=617
input=388 - output=620
input=389 - output=622
input=390 - output=625
input=391 - output=628
input=392 - output=630
input=393 - output=633
input=394 - output=635
input=395 - output=638
input=396 - output=641
input=397 - output=643
input=398 - output=646
input=399 - output=648
input=400 - output=651
input=401 - output=654
input=402 - output=656
input=403 - output=659
input=404 - output=661
input=405 - output=664
input=406 - output=667
input=407 - output=669
input=408 - output=672
input=409 - output=674
input=410 - output=677
input=411 - output=680
input=412 - output=682
input=413 - output=685
input=414 - output=687
input=415 - output=690
input=416 - output=693
input=417 - output=695
input=418 - output=698
input=419 - output=700
input=420 - output=703
input=421 - output=706
input=422 - output=708
input=423 - output=711
input=424 - output=714
input=425 - output=716
input=426 - output=719
input=427 - output=721
input=428 - output=724
input=429 - output=727
input=430 - output=729
input=431 - output=732
input=432 - output=734
input=433 - output=737
input=434 - output=740
input=435 - output=742
input=436 - output=745
input=437 - output=747
input=438 - output=750
input=439 - output=753
input=440 - output=755
input=441 - output=758
input=442 - output=760
input=443 - output=763
input=444 - output=766
input=445 - output=768
input=446 - output=771
input=447 - output=773
input=448 - output=776
input=449 - output=779
input=450 - output=781
input=451 - output=784
input=452 - output=786
input=453 - output=789
input=454 - output=792
input=455 - output=794
input=456 - output=797
input=457 - output=800
input=458 - output=800
input=459 - output=800
input=460 - output=800
input=461 - output=801
input=462 - output=801
input=463 - output=801
input=464 - output=801
input=465 - output=802
input=466 - output=802
input=467 - output=802
input=468 - output=802
input=469 - output=803
input=470 - output=803
input=471 - output=803
input=472 - output=803
input=473 - output=804
input=474 - output=804
input=475 - output=804
input=476 - output=804
input=477 - output=805
input=478 - output=805
input=479 - output=805
input=480 - output=805
input=481 - output=806
input=482 - output=806
input=483 - output=806
input=484 - output=806
input=485 - output=807
input=486 - output=807
input=487 - output=807
input=488 - output=807
input=489 - output=808
input=490 - output=808
input=491 - output=808
input=492 - output=808
input=493 - output=809
input=494 - output=809
input=495 - output=809
input=496 - output=810
input=497 - output=810
input=498 - output=810
input=499 - output=810
input=500 - output=811
input=501 - output=811
input=502 - output=811
input=503 - output=811
input=504 - output=812
input=505 - output=812
input=506 - output=812

And here is the video what my PIC outputs, the top chart is the raw data from ADC, the bottom one is with the calculated/interpolated data: https://youtu.be/Dxc60d8wDHw

Any help is highly appreciated! Thank you!

Richard Závodný
  • 11
  • 2
  • What is it that you expect from this: `uint16_t *buf = usb_get_in_buffer(1);` ? Curious, does your environment have a debugger? Does your compiler have settings to set max warnings? – ryyker Apr 12 '22 at 12:12
  • Is your 3300 figure correct? Shouldn't it be roughly 3700 to 3800? – cup Apr 25 '22 at 08:12

3 Answers3

1

You are a "victim" of arithmetic overflow and integer promotion, respectively, and that integers have different sizes on the PIC and on the PC.

All integer calculation is done with int because of integer promotion, if the operands fit in integers.

Just an example with raw = 2488:

  • a is 1988
  • b is 1988
  • c is 2700
  • d is 3500

(raw - a) gives 500 and (d - b) gives 1512.

The product of (raw - a) by (d - b) should be 756000.

On the PC you get this result, because its int uses 32 bits and has therefore a range of -2147483648 to +2147483647. After the division and the following calculations the results finally fits into uint16_t again.

But your PIC uses only 16 bits for int and unsigned int, the latter being the base type of the alias type uint16_t. Since uint16_t does not fit into an int, no integer promotion is done. However, the result overflows the range of 0 to 65535, because the product is cut down to 16 bits. The calculated result is 35104, and gives you a wrong interpolation.

The solution is to use an intermediate result with a type of int32_t or uint32_t that provides a bigger range.

the busybee
  • 10,755
  • 3
  • 13
  • 30
  • Not *all* arithmetic is done as type `int`, and on the PIC probably the OP's is not. Supposing, as seems likely, that the PIC's `uint16_t` is an alias for its `unsigned int`, there is no promotion involved there, and the result is also an `unsigned int`. But the key point here is what you say at the end: the computation overflows on the PIC, but not on Windows. – John Bollinger Apr 12 '22 at 12:34
1

In function interpolate(), this multiplication:

(raw - a) * (d - b)

produces a result of type int or unsigned int, depending on the width of those types. For an 8-bit arch, they are probably 16 bits wide (the minimum), and the result will then have type unsigned int. Again assuming 16-bit width, that result will overflow for some of your input values.

On Windows, you will find that int and unsigned int are both 32 bits wide. The result of the same multiplication will have type int (though it doesn't make a difference in this case) and it will not overflow.

There are a couple of ways you could handle this. Supposing that you want to produce the same numeric results, you can cast one of the operands to unsigned long:

    return b + (unsigned long)(raw - a) * (d - b) / (c - a);

That will probably come at a small computational cost.

Alternatively, if you were willing to accept a loss of precision then you could perform the division before the multiplication, which I would write this way for improved clarity:

    return b + (d - b) * ((raw - a) / (c - a));
John Bollinger
  • 160,171
  • 8
  • 81
  • 157
  • 1
    Instead of `unsigned long`, could use `uint32_t` as we are already in fixed-width type land. – chux - Reinstate Monica Apr 12 '22 at 13:45
  • Yes, @chux-ReinstateMonica, one could, but I see no particular advantage in doing so. We don't require a type that is exactly 32 bits wide, nor represented in two's complement form. We don't even need an unsigned type, for that matter, though I have suggested one for signedness consistency. A conforming [`unsigned`] `long` is the narrowest built-in type guaranteed to be able to represent the needed values. It's even possible for the implementation to have a conforming `long`, but not a `uint32_t`, though I doubt that's the case here. – John Bollinger Apr 12 '22 at 14:17
  • `signed long` likely will not encode `0xFFFF * 0xFFFF`. Better to use `unsigned long` than `long` for a full range solution. – chux - Reinstate Monica Apr 12 '22 at 14:34
0

As the others have analyzed, you have an arithmetic overflow problem. One solution is to use 32-bit arithmetic. Another is to keep to 16-bit arithmetic but go right down to basics.

What is 3 * 5? It is 5 + 5 + 5.

What is 10 / 3? It is how many times 3 can be subtracted from 10 until the remainder is less than 3.

10 - 3 = 7  => 1 r 7
7 - 3 = 4   => 2 r 4
4 - 3 = 1   => 3 r 1

With that in mind, you can interleave the operations. This is basically how we used to do it in assembler without multiply or divide instructions or when the differences were not nice powers of 2. It is a lot more code than a simple multiply and divide but you do not need to resort to 32-bit arithmetic. This should work on 12-bit ADCs.

uint16_t interp(uint16_t xlo, uint16_t ylo, uint16_t xhi, uint16_t yhi, uint16_t x)
{
    // yhi + (x - xlo) * (yhi - ylo) / (xhi - xlo)
    uint16_t xdiff = xhi - xlo;
    uint16_t ydiff = yhi - ylo;

    uint16_t times = x - xlo;

    // Do 16 bit interpolation without going into 32-bit arithmetic
    uint16_t remain = 0;
    uint16_t quotient = 0;
    for (int tt = 0; tt < times; ++tt)
    {
        // The 'multiply' - assume that the addition will not 
        // take the result over 0xFFFF
        remain += ydiff;

        // The 'divide'
        while (remain > xdiff)
        {
            remain -= xdiff;
            quotient += 1;
        }
    }

    // If we need to round up
    // uint16_t yhalf = ydiff >> 1;
    // if (remain > yhalf) quotient += 1;
    return ylo + quotient;
}
cup
  • 7,589
  • 4
  • 19
  • 42