comedi_calibrate/comedi_calibrate.c

   1 /*
   2    A little auto-calibration utility, for boards
   3    that support it.
   4
   5    Right now, it only supports NI E series boards,
   6    but it should be easily portable.
   7
   8    A few things need improvement here:
   9     - current system gets "close", but doesn't
  10       do any fine-tuning
  11     - no pre/post gain discrimination for the
  12       A/D zero offset.
  13     - should read (and use) the actual reference
  14       voltage value from eeprom
  15     - statistics would be nice, to show how good
  16       the calibration is.
  17     - doesn't check unipolar ranges
  18     - "alternate calibrations" would be cool--to
  19       accurately measure 0 in a unipolar range
  20     - more portable
  21  */
  22
  23 #include <stdio.h>
  24 #include <comedilib.h>
  25 #include <fcntl.h>
  26 #include <unistd.h>
  27 #include <errno.h>
  28 #include <getopt.h>
  29 #include <ctype.h>
  30 #include <math.h>
  31 #include <stdlib.h>
  32 #include <string.h>
  33
  34
  35 #define N_CALDACS 16
  36
  37 typedef struct{
  38         int subdev;
  39         int chan;
  40
  41         int maxdata;
  42         int current;
  43
  44         int type;
  45         double gain;
  46 }caldac;
  47
  48 static caldac caldacs[N_CALDACS];
  49 static int n_caldacs;
  50 static comedi_t *dev;
  51
  52 static int ad_subdev;
  53 static int eeprom_subdev;
  54
  55 double read_chan(int adc,int range);
  56 void write_caldac(comedi_t *dev,int subdev,int addr,int val);
  57 void check_gain(int ad_chan,int range);
  58 double check_gain_chan(int ad_chan,int range,int cdac);
  59
  60
  61 void update_caldac(int i);
  62 void reset_caldacs(void);
  63 void setup_caldacs(void);
  64 void cal_ni_mio_E(void);
  65 void ni_mio_ai_postgain_cal(void);
  66 void channel_dependence(int adc,int range);
  67 void caldac_dependence(int caldac);
  68 void dump_curve(int adc,int caldac);
  69 void chan_cal(int adc,int caldac,int range,double target);
  70 int read_eeprom(int addr);
  71
  72 typedef struct {
  73         int n;
  74
  75         double *y_data;
  76         double *yerr_data;
  77         double *x_data;
  78
  79         double x0;
  80         double dx;
  81         double yerr;
  82
  83         /* stats */
  84         double s1,sx,sy,sxy,sxx;
  85
  86         /* results */
  87         double ave_x;
  88         double ave_y;
  89         double slope;
  90         double err_slope;
  91         double err_ave_y;
  92         double S_min;
  93         double dof;
  94
  95 }linear_fit_t;
  96 int linear_fit_monotonic(linear_fit_t *l);
  97 double linear_fit_func_y(linear_fit_t *l,double x);
  98 double check_gain_chan_x(linear_fit_t *l,int ad_chan,int range,int cdac);
  99
 100 typedef struct{
 101         comedi_t *dev;
 102         comedi_trig t;
 103
 104         unsigned int chanlist[1];
 105
 106         int maxdata;
 107         int order;
 108         comedi_range *rng;
 109
 110         int n;
 111         double average;
 112         double stddev;
 113         double error;
 114 }new_sv_t;
 115
 116 int new_sv_measure(new_sv_t *sv);
 117 int new_sv_init(new_sv_t *sv,comedi_t *dev,int subdev,int chan,int range,int aref);
 118
 119 int main(int argc, char *argv[])
 120 {
 121         char *fn = NULL;
 122         int c;
 123         char *drivername;
 124
 125
 126         fn = "/dev/comedi0";
 127         while (1) {
 128                 c = getopt(argc, argv, "rf");
 129                 if (c == -1)
 130                         break;
 131                 switch (c) {
 132                 case 'f':
 133                         fn = argv[optind];
 134                         break;
 135                 default:
 136                         printf("bad option\n");
 137                         exit(1);
 138                 }
 139         }
 140
 141         dev = comedi_open(fn);
 142         if (dev == NULL ) {
 143                 perror(fn);
 144                 exit(0);
 145         }
 146
 147         ad_subdev=0;
 148         setup_caldacs();
 149
 150         eeprom_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_MEMORY,0);
 151
 152         drivername=comedi_get_driver_name(dev);
 153
 154         if(!strcmp(drivername,"atmio-E") || !strcmp(drivername,"pcimio-E"))
 155                 cal_ni_mio_E();
 156
 157         return 0;
 158 }
 159
 160
 161 void cal_ni_mio_E(void)
 162 {
 163         char *boardname;
 164         double ref;
 165         int i;
 166
 167         boardname=comedi_get_board_name(dev);
 168
 169         if(!strcmp(boardname,"at-mio-16e-1") ||
 170            !strcmp(boardname,"at-mio-16e-2") ||
 171            !strcmp(boardname,"at-mio-64e-3") ||
 172            !strcmp(boardname,"pci-mio-16e-1")){
 173 /*
 174    layout
 175
 176    0    AI pre-gain offset      1.5e-6
 177    1    AI post-gain offset     8.1e-4
 178    2    AI unipolar offset      7.9e-4
 179    3    AI gain
 180    4    AO
 181    5    AO
 182    6    AO
 183    7    A0
 184    8    AO
 185    9    AO
 186    10   analog trigger
 187    11   unknown
 188  */
 189                 printf("last factory calibration %02d/%02d/%02d\n",
 190                         read_eeprom(508),read_eeprom(507),read_eeprom(506));
 191
 192                 printf("lsb=%d msb=%d\n",read_eeprom(425),read_eeprom(426));
 193
 194                 ref=5.000+(1e-6*(read_eeprom(425)+read_eeprom(426)));
 195                 printf("ref=%g\n",ref);
 196
 197                 reset_caldacs();
 198
 199                 printf("postgain offset\n");
 200                 ni_mio_ai_postgain_cal();
 201
 202                 printf("pregain offset\n");
 203                 chan_cal(0,0,7,0.0);
 204                 chan_cal(0,0,7,0.0);
 205
 206                 printf("unipolar offset\n");
 207                 chan_cal(0,2,8,0.0);
 208                 chan_cal(0,2,8,0.0);
 209
 210                 printf("gain offset\n");
 211                 chan_cal(5,3,0,5.0);
 212                 chan_cal(5,3,0,5.0);
 213
 214                 return;
 215         }
 216         if(!strcmp(boardname,"at-mio-16e-10")){
 217 /*
 218    layout
 219
 220    0    AI pre-gain offset      1.5e-6
 221    1    AI post-gain offset     8.1e-4
 222    2    AI unipolar offset      7.9e-4
 223    3    AI gain                 3.5e-4
 224    4    AO
 225    5    AO
 226    6    AO
 227    7    AO
 228    8    AO
 229    9    AO
 230    10   AI pre-gain offset      6.4e-5
 231    11   unknown
 232  */
 233                 printf("last factory calibration %02d/%02d/%02d\n",
 234                         read_eeprom(508),read_eeprom(507),read_eeprom(506));
 235
 236                 printf("lsb=%d msb=%d\n",read_eeprom(423),read_eeprom(424));
 237
 238                 ref=5.000+(0.001*(read_eeprom(423)+read_eeprom(424)));
 239                 printf("ref=%g\n",ref);
 240
 241                 reset_caldacs();
 242
 243                 printf("postgain offset\n");
 244                 ni_mio_ai_postgain_cal();
 245
 246                 printf("pregain offset\n");
 247                 chan_cal(0,10,7,0.0);
 248                 chan_cal(0,0,7,0.0);
 249                 chan_cal(0,0,7,0.0);
 250
 251                 printf("unipolar offset\n");
 252                 chan_cal(0,2,8,0.0);
 253                 chan_cal(0,2,8,0.0);
 254
 255                 printf("gain offset\n");
 256                 chan_cal(5,3,0,5.0);
 257                 chan_cal(5,3,0,5.0);
 258
 259                 printf("results (offset)\n");
 260                 for(i=0;i<16;i++){
 261                         read_chan(0,i);
 262                 }
 263
 264                 return;
 265         }
 266
 267         {
 268                 int n_ranges;
 269
 270         reset_caldacs();
 271         printf("please send this output to <ds@stm.lbl.gov>\n");
 272         printf("%s\n",comedi_get_board_name(dev));
 273
 274         n_ranges=comedi_get_n_ranges(dev,ad_subdev,0);
 275
 276         printf("channel dependence 0 range 0\n");
 277         channel_dependence(0,0);
 278
 279         printf("channel dependence 0 range %d\n",n_ranges/2-1);
 280         channel_dependence(0,n_ranges/2-1);
 281
 282         printf("channel dependence 0 range %d\n",n_ranges/2);
 283         channel_dependence(0,n_ranges/2);
 284
 285         printf("channel dependence 5 range 0\n");
 286         channel_dependence(5,0);
 287         }
 288 }
 289
 290
 291 void ni_mio_ai_postgain_cal(void)
 292 {
 293         linear_fit_t l;
 294         double offset_r0;
 295         double offset_r7;
 296         double gain;
 297         double a;
 298
 299         check_gain_chan_x(&l,0,0,1);
 300         offset_r0=linear_fit_func_y(&l,caldacs[1].current);
 301         printf("offset r0 %g\n",offset_r0);
 302
 303         check_gain_chan_x(&l,0,7,1);
 304         offset_r7=linear_fit_func_y(&l,caldacs[1].current);
 305         printf("offset r7 %g\n",offset_r7);
 306
 307         gain=l.slope;
 308
 309         a=(offset_r0-offset_r7)/(200.0-1.0);
 310         a=caldacs[1].current-a/gain;
 311
 312         printf("%g\n",a);
 313
 314         caldacs[1].current=rint(a);
 315         update_caldac(1);
 316 }
 317
 318 void chan_cal(int adc,int cdac,int range,double target)
 319 {
 320         linear_fit_t l;
 321         double offset;
 322         double gain;
 323         double a;
 324
 325         check_gain_chan_x(&l,adc,range,cdac);
 326         offset=linear_fit_func_y(&l,caldacs[cdac].current);
 327         printf("offset %g\n",offset);
 328         gain=l.slope;
 329
 330         a=caldacs[cdac].current+(target-offset)/gain;
 331         printf("%g\n",a);
 332
 333         caldacs[cdac].current=rint(a);
 334         update_caldac(cdac);
 335 }
 336
 337
 338
 339 void channel_dependence(int adc,int range)
 340 {
 341         int i;
 342         double gain;
 343
 344         for(i=0;i<n_caldacs;i++){
 345                 gain=check_gain_chan(adc,range,i);
 346         }
 347 }
 348
 349 void caldac_dependence(int caldac)
 350 {
 351         int i;
 352         double gain;
 353
 354         for(i=0;i<8;i++){
 355                 gain=check_gain_chan(i,0,caldac);
 356         }
 357 }
 358
 359 int dump_flag;
 360
 361 void dump_curve(int adc,int caldac)
 362 {
 363         linear_fit_t l;
 364
 365         dump_flag=1;
 366         check_gain_chan_x(&l,adc,0,caldac);
 367         dump_flag=0;
 368 }
 369
 370
 371 void setup_caldacs(void)
 372 {
 373         int s,n_chan,i;
 374
 375         s=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_CALIB,0);
 376         if(s<0){
 377                 printf("no calibration subdevice\n");
 378                 return;
 379         }
 380         //printf("calibration subdevice is %d\n",s);
 381
 382         n_chan=comedi_get_n_channels(dev,s);
 383
 384         for(i=0;i<n_chan;i++){
 385                 caldacs[i].subdev=s;
 386                 caldacs[i].chan=i;
 387                 caldacs[i].maxdata=comedi_get_maxdata(dev,s,i);
 388 /* XXX */
 389 caldacs[i].maxdata=255;
 390                 caldacs[i].current=0;
 391                 //printf("caldac %d, %d\n",i,caldacs[i].maxdata);
 392         }
 393
 394         n_caldacs=n_chan;
 395 }
 396
 397 void reset_caldacs(void)
 398 {
 399         int i;
 400
 401         for(i=0;i<n_caldacs;i++){
 402                 caldacs[i].current=caldacs[i].maxdata/2;
 403                 update_caldac(i);
 404         }
 405 }
 406
 407 void update_caldac(int i)
 408 {
 409         int ret;
 410
 411         //printf("update %d %d %d\n",caldacs[i].subdev,caldacs[i].chan,caldacs[i].current);
 412         //ret=comedi_data_write(dev,caldacs[i].subdev,caldacs[i].chan,0,0,caldacs[i].current);
 413         write_caldac(dev,caldacs[i].subdev,caldacs[i].chan,caldacs[i].current);
 414         ret=0;
 415         if(ret<0)
 416                 printf("comedi_data_write() error\n");
 417 }
 418
 419
 420 void check_gain(int ad_chan,int range)
 421 {
 422         int i;
 423
 424         for(i=0;i<n_caldacs;i++){
 425                 check_gain_chan(ad_chan,range,i);
 426         }
 427 }
 428
 429 double check_gain_chan(int ad_chan,int range,int cdac)
 430 {
 431         linear_fit_t l;
 432
 433         return check_gain_chan_x(&l,ad_chan,range,cdac);
 434 }
 435
 436 double check_gain_chan_x(linear_fit_t *l,int ad_chan,int range,int cdac)
 437 {
 438         int orig,i,n;
 439         new_sv_t sv;
 440         double sum_err;
 441         int sum_err_count=0;
 442
 443         n=caldacs[cdac].maxdata+1;
 444         memset(l,0,sizeof(*l));
 445         l->y_data=malloc(n*sizeof(double));
 446
 447         orig=caldacs[cdac].current;
 448
 449         new_sv_init(&sv,dev,0,ad_chan,range,AREF_OTHER);
 450
 451         sum_err=0;
 452         for(i=0;i<n;i++){
 453                 caldacs[cdac].current=i;
 454                 update_caldac(cdac);
 455
 456                 new_sv_measure(&sv);
 457
 458                 l->y_data[i]=sv.average;
 459                 if(!isnan(sv.average)){
 460                         sum_err+=sv.error;
 461                         sum_err_count++;
 462                 }
 463         }
 464
 465         caldacs[cdac].current=orig;
 466         update_caldac(cdac);
 467
 468         l->yerr=sum_err/sum_err_count;
 469         l->n=n;
 470         l->dx=1;
 471         l->x0=0;
 472
 473         linear_fit_monotonic(l);
 474
 475         printf("caldac[%d] gain=%g V/bit err=%g S_min=%g dof=%g\n",
 476                 cdac,l->slope,l->err_slope,l->S_min,l->dof);
 477         //printf("--> %g\n",fabs(l.slope/l.err_slope));
 478
 479         if(dump_flag){
 480                 for(i=0;i<n;i++){
 481                         printf("%d %g\n",i,l->y_data[i]);
 482                 }
 483         }
 484
 485
 486         free(l->y_data);
 487
 488         return l->slope;
 489 }
 490
 491
 492
 493
 494 /* helpers */
 495
 496
 497 int read_eeprom(int addr)
 498 {
 499         unsigned int data=0;
 500
 501         comedi_data_read(dev,eeprom_subdev,addr,0,0,&data);
 502
 503         return data;
 504 }
 505
 506 double read_chan(int adc,int range)
 507 {
 508         int n;
 509         new_sv_t sv;
 510
 511         new_sv_init(&sv,dev,0,adc,range,AREF_OTHER);
 512         sv.order=10;
 513         n=new_sv_measure(&sv);
 514
 515         printf("chan=%d ave=%g error=%g\n",adc,sv.average,sv.error);
 516
 517         return sv.average;
 518 }
 519
 520 void write_caldac(comedi_t *dev,int subdev,int addr,int val)
 521 {
 522         comedi_trig it;
 523         unsigned int chan=CR_PACK(addr,0,0);
 524         sampl_t data=val;
 525
 526         memset(&it,0,sizeof(it));
 527
 528         it.subdev = subdev;
 529         it.n_chan = 1;
 530         it.chanlist = &chan;
 531         it.data = &data;
 532         it.n = 1;
 533
 534         if(comedi_trigger(dev,&it)<0){
 535                 perror("write_caldac");
 536                 exit(0);
 537         }
 538 }
 539
 540
 541
 542 int new_sv_init(new_sv_t *sv,comedi_t *dev,int subdev,int chan,int range,int aref)
 543 {
 544         memset(sv,0,sizeof(*sv));
 545
 546         sv->t.subdev=subdev;
 547         sv->t.flags=TRIG_DITHER;
 548         sv->t.n_chan=1;
 549         sv->t.chanlist=sv->chanlist;
 550
 551         sv->chanlist[0]=CR_PACK(chan,range,aref);
 552
 553         sv->maxdata=comedi_get_maxdata(dev,subdev,chan);
 554         sv->rng=comedi_get_range(dev,subdev,chan,range);
 555
 556         sv->order=7;
 557
 558         return 0;
 559 }
 560
 561 int new_sv_measure(new_sv_t *sv)
 562 {
 563         sampl_t *data;
 564         int n,i,ret;
 565
 566         double a,x,s,s2;
 567
 568         n=1<<sv->order;
 569
 570         data=malloc(sizeof(sampl_t)*n);
 571
 572         for(i=0;i<n;){
 573                 sv->t.data=data+i;
 574                 sv->t.n=n-i;
 575                 ret=comedi_trigger(dev,&sv->t);
 576                 if(ret<0)goto out;
 577                 i+=ret;
 578         }
 579
 580         s=0;
 581         s2=0;
 582         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
 583         for(i=0;i<n;i++){
 584                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
 585                 s+=x-a;
 586                 s2+=(x-a)*(x-a);
 587         }
 588         sv->average=a+s/n;
 589         sv->stddev=sqrt(n*s2-s*s)/n;
 590         sv->error=sv->stddev/sqrt(n);
 591
 592         ret=n;
 593
 594 out:
 595         free(data);
 596
 597         return ret;
 598 }
 599
 600 int new_sv_measure_order(new_sv_t *sv,int order)
 601 {
 602         sampl_t *data;
 603         int n,i,ret;
 604         double a,x,s,s2;
 605
 606         n=1<<order;
 607
 608         data=malloc(sizeof(sampl_t)*n);
 609
 610         for(i=0;i<n;){
 611                 sv->t.data=data+i;
 612                 sv->t.n=n-i;
 613                 ret=comedi_trigger(dev,&sv->t);
 614                 if(ret<0)goto out;
 615                 i+=ret;
 616         }
 617
 618         s=0;
 619         s2=0;
 620         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
 621         for(i=0;i<n;i++){
 622                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
 623                 s+=x-a;
 624                 s2+=(x-a)*(x-a);
 625         }
 626         sv->average=a+s/n;
 627         sv->stddev=sqrt(n*s2-s*s)/n;
 628         sv->error=sv->stddev/sqrt(n);
 629
 630         ret=n;
 631
 632 out:
 633         free(data);
 634
 635         return ret;
 636 }
 637
 638
 639
 640 /* linear fitting */
 641
 642 int calculate_residuals(linear_fit_t *l);
 643
 644 int linear_fit_monotonic(linear_fit_t *l)
 645 {
 646         double x,y;
 647         double sxp;
 648         int i;
 649
 650         l->s1=0;
 651         l->sx=0;
 652         l->sy=0;
 653         l->sxy=0;
 654         l->sxx=0;
 655         for(i=0;i<l->n;i++){
 656                 x=l->x0+i*l->dx;
 657                 y=l->y_data[i];
 658
 659                 if(isnan(y))continue;
 660
 661                 l->s1+=1;
 662                 l->sx+=x;
 663                 l->sy+=y;
 664                 l->sxy+=x*y;
 665                 l->sxx+=x*x;
 666         }
 667         sxp=l->sxx-l->sx*l->sx/l->s1;
 668
 669         l->ave_x=l->sx/l->s1;
 670         l->ave_y=l->sy/l->s1;
 671         l->slope=(l->s1*l->sxy-l->sx*l->sy)/(l->s1*l->sxx-l->sx*l->sx);
 672         l->err_slope=l->yerr/sqrt(sxp);
 673         l->err_ave_y=l->yerr/sqrt(l->s1);
 674
 675         calculate_residuals(l);
 676
 677         return 0;
 678 }
 679
 680 int calculate_residuals(linear_fit_t *l)
 681 {
 682         double x,y;
 683         double res,sum_res2;
 684         int i;
 685
 686         sum_res2=0;
 687         for(i=0;i<l->n;i++){
 688                 x=l->x0+i*l->dx-l->ave_x;
 689                 y=l->y_data[i];
 690
 691                 if(isnan(y))continue;
 692
 693                 res=l->ave_y+l->slope*x-y;
 694
 695                 sum_res2+=res*res;
 696         }
 697         l->S_min=sum_res2/(l->yerr*l->yerr);
 698         l->dof=l->s1-2;
 699
 700         return 0;
 701 }
 702
 703 double linear_fit_func_y(linear_fit_t *l,double x)
 704 {
 705         return l->ave_y+l->slope*(x-l->ave_x);
 706 }
 707