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 #define _GNU_SOURCE
  24
  25 #include <stdio.h>
  26 #include <comedilib.h>
  27 #include <fcntl.h>
  28 #include <unistd.h>
  29 #include <errno.h>
  30 #include <getopt.h>
  31 #include <ctype.h>
  32 #include <math.h>
  33 #include <stdlib.h>
  34 #include <string.h>
  35
  36
  37 #define N_CALDACS 32
  38
  39 typedef struct{
  40         int subdev;
  41         int chan;
  42
  43         int maxdata;
  44         int current;
  45
  46         int type;
  47         double gain;
  48 }caldac;
  49
  50 static caldac caldacs[N_CALDACS];
  51 static int n_caldacs;
  52 static comedi_t *dev;
  53
  54 static int ad_subdev;
  55 static int eeprom_subdev;
  56
  57 double read_chan(int adc,int range);
  58 void write_caldac(comedi_t *dev,int subdev,int addr,int val);
  59 void check_gain(int ad_chan,int range);
  60 double check_gain_chan(int ad_chan,int range,int cdac);
  61
  62 int verbose = 1;
  63
  64
  65 void update_caldac(int i);
  66 void reset_caldacs(void);
  67 void setup_caldacs(void);
  68 void cal_ni_mio_E(void);
  69 void ni_mio_ai_postgain_cal(void);
  70 void ni_mio_ai_postgain_cal_2(int chan,int dac,int range_lo,int range_hi,double gain);
  71 void channel_dependence(int adc,int range);
  72 void caldac_dependence(int caldac);
  73 void dump_curve(int adc,int caldac);
  74 void chan_cal(int adc,int caldac,int range,double target);
  75 int read_eeprom(int addr);
  76
  77 typedef struct {
  78         int n;
  79
  80         double *y_data;
  81         double *yerr_data;
  82         double *x_data;
  83
  84         double x0;
  85         double dx;
  86         double yerr;
  87
  88         /* stats */
  89         double s1,sx,sy,sxy,sxx;
  90
  91         double min,max;
  92
  93         /* results */
  94         double ave_x;
  95         double ave_y;
  96         double slope;
  97         double err_slope;
  98         double err_ave_y;
  99         double S_min;
 100         double dof;
 101
 102 }linear_fit_t;
 103 int linear_fit_monotonic(linear_fit_t *l);
 104 double linear_fit_func_y(linear_fit_t *l,double x);
 105 double check_gain_chan_x(linear_fit_t *l,int ad_chan,int range,int cdac);
 106
 107 typedef struct{
 108         comedi_t *dev;
 109
 110         int maxdata;
 111         int order;
 112         int aref;
 113         int range;
 114         int subd;
 115         int chan;
 116
 117         comedi_range *rng;
 118
 119         int n;
 120         double average;
 121         double stddev;
 122         double error;
 123 }new_sv_t;
 124
 125 int new_sv_measure(new_sv_t *sv);
 126 int new_sv_init(new_sv_t *sv,comedi_t *dev,int subdev,int chan,int range,int aref);
 127
 128 struct board_struct{
 129         char *name;
 130         void (*calibrate)(void);
 131 };
 132
 133 void cal_ni_16e_1(void);
 134 void cal_ni_16e_10(void);
 135 void cal_ni_16xe_50(void);
 136 void cal_ni_16xe_10(void);
 137 void cal_ni_6023e(void);
 138 void cal_ni_daqcard_ai_16xe_50(void);
 139 void cal_ni_unknown(void);
 140
 141 struct board_struct boards[]={
 142         { "at-mio-16e-1",       cal_ni_16e_1 },
 143         { "at-mio-16e-2",       cal_ni_16e_1 },
 144         { "at-mio-16e-10",      cal_ni_16e_10 },
 145         { "at-mio-16de-10",     cal_ni_unknown },
 146         { "at-mio-64e-3",       cal_ni_16e_1 },
 147         { "at-mio-16xe-50",     cal_ni_unknown },
 148         { "at-mio-16xe-10",     cal_ni_unknown },
 149         { "at-ai-16xe-10",      cal_ni_unknown },
 150         { "pci-mio-16xe-50",    cal_ni_16xe_50 },
 151         { "pci-mio-16xe-10",    cal_ni_16xe_10 },
 152         { "pxi-6030e",          cal_ni_unknown },
 153         { "pci-mio-16e-1",      cal_ni_16e_1 },
 154         { "pci-mio-16e-4",      cal_ni_unknown },
 155         { "pxi-6040e",          cal_ni_unknown },
 156         { "pci-6031e",          cal_ni_unknown },
 157         { "pci-6032e",          cal_ni_unknown },
 158         { "pci-6033e",          cal_ni_unknown },
 159         { "pci-6071e",          cal_ni_unknown },
 160         { "pci-6023e",          cal_ni_6023e },
 161         { "pci-6024e",          cal_ni_unknown },
 162         { "pci-6025e",          cal_ni_unknown },
 163         { "pxi-6025e",          cal_ni_unknown },
 164         { "pci-6034e",          cal_ni_unknown },
 165         { "pci-6035e",          cal_ni_unknown },
 166         { "pci-6052e",          cal_ni_unknown },
 167         { "pci-6110e",          cal_ni_unknown },
 168         { "pci-6111e",          cal_ni_unknown },
 169 //      { "pci-6711",           cal_ni_unknown },
 170 //      { "pci-6713",           cal_ni_unknown },
 171         { "pxi-6071e",          cal_ni_unknown },
 172         { "pxi-6070e",          cal_ni_unknown },
 173         { "pxi-6052e",          cal_ni_unknown },
 174         { "DAQCard-ai-16xe-50", cal_ni_daqcard_ai_16xe_50 },
 175         { "DAQCard-ai-16e-4",   cal_ni_unknown },
 176         { "DAQCard-6062e",      cal_ni_unknown },
 177         { "DAQCard-6024e",      cal_ni_unknown },
 178 };
 179 #define n_boards (sizeof(boards)/sizeof(boards[0]))
 180
 181 int main(int argc, char *argv[])
 182 {
 183         char *fn = NULL;
 184         int c;
 185         char *drivername;
 186         char *devicename;
 187         int i;
 188
 189         fn = "/dev/comedi0";
 190         while (1) {
 191                 c = getopt(argc, argv, "f:vq");
 192                 if (c == -1)
 193                         break;
 194                 switch (c) {
 195                 case 'f':
 196                         fn = optarg;
 197                         break;
 198                 case 'v':
 199                         verbose++;
 200                         break;
 201                 case 'q':
 202                         verbose--;
 203                         break;
 204                 default:
 205                         printf("bad option\n");
 206                         exit(1);
 207                 }
 208         }
 209
 210         dev = comedi_open(fn);
 211         if (dev == NULL ) {
 212                 perror(fn);
 213                 exit(0);
 214         }
 215
 216         ad_subdev=0;
 217         setup_caldacs();
 218
 219         eeprom_subdev=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_MEMORY,0);
 220
 221         drivername=comedi_get_driver_name(dev);
 222         devicename=comedi_get_board_name(dev);
 223
 224         for(i=0;i<n_boards;i++){
 225                 if(!strcmp(boards[i].name,devicename)){
 226                         boards[i].calibrate();
 227                         break;
 228                 }
 229         }
 230
 231         printf("device %s unknown\n",devicename);
 232
 233         return 0;
 234 }
 235
 236 double ni_get_reference(int lsb_loc,int msb_loc)
 237 {
 238         int lsb,msb;
 239         int uv;
 240         double ref;
 241
 242         lsb=read_eeprom(lsb_loc);
 243         msb=read_eeprom(msb_loc);
 244         printf("lsb=%d msb=%d\n",read_eeprom(425),read_eeprom(426));
 245
 246         uv=lsb | (msb<<8);
 247         if(uv>=0x8000)uv-=0x10000;
 248         ref=5.000+1.0e-6*uv;
 249         printf("ref=%g\n",ref);
 250
 251         return ref;
 252 }
 253
 254 void cal_ni_16e_1(void)
 255 {
 256         double ref;
 257
 258         reset_caldacs();
 259
 260 /*
 261    layout
 262
 263    0    AI pre-gain offset      1.5e-6
 264    1    AI post-gain offset     8.1e-4
 265    2    AI unipolar offset      7.9e-4
 266    3    AI gain
 267    4    AO 0    -1.2e-4         -1.2e-4
 268    5    AO 0    -8.0e-4         -8.0e-4
 269    6    AO 0    1.9e-4          -3.8e-7
 270    7    AO 1    -8.0e-5         -1.2e-4
 271    8    AO 1    -7.9e-4         -7.9e-4
 272    9    AO 1    1.9e-4          3.0e-7
 273    10   analog trigger
 274    11   unknown
 275  */
 276         printf("last factory calibration %02d/%02d/%02d\n",
 277                 read_eeprom(508),read_eeprom(507),read_eeprom(506));
 278
 279         ref=ni_get_reference(425,426);
 280
 281         printf("postgain offset\n");
 282         ni_mio_ai_postgain_cal();
 283
 284         printf("pregain offset\n");
 285         chan_cal(0,0,7,0.0);
 286         chan_cal(0,0,7,0.0);
 287
 288         printf("unipolar offset\n");
 289         chan_cal(0,2,8,0.0);
 290         chan_cal(0,2,8,0.0);
 291
 292         printf("gain offset\n");
 293         chan_cal(5,3,0,5.0);
 294         chan_cal(5,3,0,5.0);
 295
 296         printf("ao 0 offset\n");
 297         comedi_data_write(dev,1,0,0,0,2048);
 298         chan_cal(2,4,0,0.0);
 299         chan_cal(2,5,0,0.0);
 300
 301         printf("ao 0 gain\n");
 302         comedi_data_write(dev,1,0,0,0,3072);
 303         chan_cal(6,6,0,0.0);
 304         chan_cal(6,6,0,0.0);
 305         comedi_data_write(dev,1,0,0,0,2048);
 306 }
 307
 308
 309 void cal_ni_16e_10(void)
 310 {
 311         double ref;
 312         int i;
 313
 314 /*
 315    layout
 316
 317    0    AI pre-gain offset      1.5e-6
 318    1    AI post-gain offset     8.1e-4
 319    2    AI unipolar offset      7.9e-4
 320    3    AI gain                 3.5e-4
 321    4    AO
 322    5    AO
 323    6    AO
 324    7    AO
 325    8    AO
 326    9    AO
 327    10   AI pre-gain offset      6.4e-5
 328    11   unknown
 329  */
 330         printf("last factory calibration %02d/%02d/%02d\n",
 331                 read_eeprom(508),read_eeprom(507),read_eeprom(506));
 332
 333         ref=ni_get_reference(423,424);
 334
 335         reset_caldacs();
 336
 337         printf("postgain offset\n");
 338         ni_mio_ai_postgain_cal();
 339
 340         printf("pregain offset\n");
 341         chan_cal(0,10,7,0.0);
 342         chan_cal(0,0,7,0.0);
 343         chan_cal(0,0,7,0.0);
 344
 345         printf("unipolar offset\n");
 346         chan_cal(0,2,8,0.0);
 347         chan_cal(0,2,8,0.0);
 348
 349         printf("gain offset\n");
 350         chan_cal(5,3,0,5.0);
 351         chan_cal(5,3,0,5.0);
 352
 353         printf("results (offset)\n");
 354         for(i=0;i<16;i++){
 355                 read_chan(i,0);
 356         }
 357 }
 358
 359 void cal_ni_16xe_10(void)
 360 {
 361         double ref;
 362         int i;
 363
 364 /*
 365  * results of channel dependency test:
 366  *
 367  *              [0]     [1]     [2]     [3]     [8]
 368  * offset, lo                   1.9e-4* 2.2e-6  2.4e-7
 369  * offset, hi                   2.0e-6* 2.1e-8  2.7e-7
 370  * offset, unip                 1.9e-4  2.1e-6  3.9e-7
 371  * ref          -2.3e-5*-1.3e-6*1.9e-4* 2.1e-6* 3.2e-7
 372  *
 373  * thus, 2,3 are postgain offset, 8 is pregain, and
 374  * 0,1 is gain.  Note the suspicious lack of unipolar
 375  * offset.
 376  *
 377  * layout
 378  *
 379  * 0    AI gain                 -2.3e-5
 380  * 1    AI gain                 -1.3e-6
 381  * 2    AI postgain offset      1.9e-4
 382  * 3    AI postgain offset      2.2e-6
 383  * 4    AO
 384  * 5    AO
 385  * 6    AO
 386  * 7    AO
 387  * 8    AI pregain offset       2.4e-7
 388  * 9    unknown
 389  * 10   unknown
 390  */
 391         printf("last factory calibration %02d/%02d/%02d\n",
 392                 read_eeprom(508),read_eeprom(507),read_eeprom(506));
 393
 394         ref=ni_get_reference(430,431);
 395
 396         reset_caldacs();
 397
 398         printf("postgain offset\n");
 399         ni_mio_ai_postgain_cal_2(0,2,0,6,100.0);
 400         ni_mio_ai_postgain_cal_2(0,3,0,6,100.0);
 401
 402         printf("pregain offset\n");
 403         chan_cal(0,8,6,0.0);
 404         chan_cal(0,8,6,0.0);
 405
 406         //printf("unipolar offset\n");
 407         //chan_cal(0,2,8,0.0);
 408         //chan_cal(0,2,8,0.0);
 409
 410         printf("gain offset\n");
 411         chan_cal(5,0,0,5.0);
 412         chan_cal(5,1,0,5.0);
 413         chan_cal(5,1,0,5.0);
 414
 415         printf("results (offset)\n");
 416         for(i=0;i<16;i++){
 417                 read_chan(i,0);
 418         }
 419 }
 420
 421 void cal_ni_daqcard_ai_16xe_50(void)
 422 {
 423         double ref;
 424         int i;
 425
 426 /*
 427  * results of channel dependency test:
 428  *
 429  *              [0]     [1]     [2]     [3]     [8]
 430  * offset, lo   -2.2e-6         1.5e-4*         2.5e-7
 431  * offset, hi                   7.8e-7*         1.3e-7
 432  * offset, unip 7.4e-4  1.1e-5  1.5e-4          5.5e-7
 433  * ref          -3.7e-4 -5.4e-6 1.5e-4*         5.5e-7
 434  *
 435  * thus, 2 is postgain offset, 8 is pregain, 0 is
 436  * unipolar offset, 1 is gain
 437  *
 438  * layout
 439  *
 440  * 0    AI unipolar offset      7.4e-4
 441  * 1    AI gain                 -5.4e-6
 442  * 2    AI postgain offset      1.5e-4
 443  * 3    unknown
 444  * 4    AO
 445  * 5    AO
 446  * 6    AO
 447  * 7    AO
 448  * 8    AI pregain offset       2.5e-7
 449  * 9    unknown
 450  * 10   unknown
 451  */
 452         printf("last factory calibration %02d/%02d/%02d\n",
 453                 read_eeprom(508),read_eeprom(507),read_eeprom(506));
 454
 455         ref=ni_get_reference(446,447);
 456
 457         reset_caldacs();
 458
 459         printf("postgain offset\n");
 460         ni_mio_ai_postgain_cal_2(0,2,0,3,100.0);
 461
 462         printf("pregain offset\n");
 463         chan_cal(0,8,3,0.0);
 464         chan_cal(0,8,3,0.0);
 465
 466         printf("unipolar offset\n");
 467         chan_cal(0,0,4,0.0);
 468         chan_cal(0,0,4,0.0);
 469
 470         printf("gain offset\n");
 471         chan_cal(5,1,0,ref);
 472         chan_cal(5,1,0,ref);
 473
 474         printf("results (offset)\n");
 475         for(i=0;i<8;i++){
 476                 read_chan(0,i);
 477         }
 478 }
 479
 480 void cal_ni_16xe_50(void)
 481 {
 482         double ref;
 483         int i;
 484
 485 /*
 486  * results of channel dependency test:
 487  *
 488  *              [0]     [1]     [2]     [3]     [8]
 489  * offset, lo                   1.6e-5          2.0e-7
 490  * offset, hi                   1.6e-7          1.8e-7
 491  * offset, unip
 492  * ref          -4.5e-5 -2.9e-6 1.6e-5*         5.5e-7
 493  *
 494  * thus, 2 is postgain offset, 8 is pregain, 0 is
 495  * unipolar offset, 1 is gain
 496  *
 497  * layout
 498  *
 499  * 0    AI unipolar offset      7.4e-4
 500  * 1    AI gain                 -5.4e-6
 501  * 2    AI postgain offset      1.5e-4
 502  * 3    unknown
 503  * 4    AO
 504  * 5    AO
 505  * 6    AO
 506  * 7    AO
 507  * 8    AI pregain offset       2.5e-7
 508  * 9    unknown
 509  * 10   unknown
 510  */
 511         printf("last factory calibration %02d/%02d/%02d\n",
 512                 read_eeprom(508),read_eeprom(507),read_eeprom(506));
 513
 514         ref=ni_get_reference(437,438);
 515
 516         reset_caldacs();
 517
 518         printf("postgain offset\n");
 519         ni_mio_ai_postgain_cal_2(0,2,0,3,100.0);
 520
 521         printf("pregain offset\n");
 522         chan_cal(0,8,3,0.0);
 523         chan_cal(0,8,3,0.0);
 524
 525 #if 0
 526         printf("unipolar offset\n");
 527         chan_cal(0,0,4,0.0);
 528         chan_cal(0,0,4,0.0);
 529 #endif
 530
 531         printf("gain offset\n");
 532         chan_cal(5,0,0,5.0);
 533         chan_cal(5,1,0,5.0);
 534         chan_cal(5,1,0,5.0);
 535
 536         printf("results (offset)\n");
 537         for(i=0;i<16;i++){
 538                 read_chan(0,i);
 539         }
 540 }
 541
 542 void cal_ni_6023e(void)
 543 {
 544         double ref;
 545         int i;
 546
 547 /*
 548  * results of channel dependency test:
 549  *
 550  *              [0]     [1]     [3]     [10]
 551  * offset, lo   -2.8e-9 -7.6e-4
 552  * offset, hi   -2.0e-6 -3.8e-6 -1.4e-6
 553  * offset, unip         1.0e-1*
 554  * ref          -7.6e-7 -7.6e-4 -5.6e-4 -6.2e-8
 555  * ref2         -6.3e-8 -7.5e-4 -5.6e-4 -1.5e-8
 556  *
 557  * 0 is pregain offset
 558  * 1 is postgain offset
 559  * 3 is gain
 560  *
 561  * layout
 562  *
 563  * 0    AI pregain offset       -2.0e-6
 564  * 1    AI postgain offset      -7.6e-4
 565  * 2    unknown
 566  * 3    AI gain                 -5.6e-4
 567  * 4    AO
 568  * 5    AO
 569  * 6    AO
 570  * 7    AO
 571  * 8    unknown
 572  * 9    unknown
 573  * 10   AI                      ?
 574  * 11   unknown
 575  */
 576         int offset_ad = 0;
 577         //int unipolar_offset_ad = 1;
 578         int gain_ad = 5;
 579         int pregain_offset_dac = 0;
 580         int postgain_offset_dac = 1;
 581         int gain_dac = 3;
 582
 583         printf("last factory calibration %02d/%02d/%02d\n",
 584                 read_eeprom(508),read_eeprom(507),read_eeprom(506));
 585
 586         ref=ni_get_reference(444,443);
 587
 588         reset_caldacs();
 589
 590         printf("postgain offset\n");
 591         ni_mio_ai_postgain_cal_2(offset_ad,postgain_offset_dac,0,3,200.0);
 592
 593         printf("pregain offset\n");
 594         chan_cal(offset_ad,pregain_offset_dac,3,0.0);
 595         chan_cal(offset_ad,pregain_offset_dac,3,0.0);
 596
 597         printf("gain offset\n");
 598         chan_cal(gain_ad,gain_dac,0,5.0);
 599         chan_cal(gain_ad,gain_dac,0,5.0);
 600
 601         printf("results (offset)\n");
 602         for(i=0;i<16;i++){
 603                 read_chan(0,i);
 604         }
 605 }
 606
 607 void cal_ni_unknown(void)
 608 {
 609         int n_ranges;
 610
 611         reset_caldacs();
 612         printf("Please send this output to <ds@schleef.org>\n");
 613         printf("$Id$\n");
 614         printf("Device name: %s\n",comedi_get_board_name(dev));
 615         printf("Comedi version: %d.%d.%d\n",
 616                 (comedi_get_version_code(dev)>>16)&0xff,
 617                 (comedi_get_version_code(dev)>>8)&0xff,
 618                 (comedi_get_version_code(dev))&0xff);
 619
 620         n_ranges=comedi_get_n_ranges(dev,ad_subdev,0);
 621
 622         /* 0 offset, low gain */
 623         printf("channel dependence 0 range 0\n");
 624         channel_dependence(0,0);
 625
 626         /* 0 offset, high gain */
 627         printf("channel dependence 0 range %d\n",n_ranges/2-1);
 628         channel_dependence(0,n_ranges/2-1);
 629
 630         /* unip/bip offset */
 631         printf("channel dependence 0 range %d\n",n_ranges/2);
 632         channel_dependence(0,n_ranges/2);
 633
 634         /* voltage reference */
 635         printf("channel dependence 5 range 0\n");
 636         channel_dependence(5,0);
 637
 638 }
 639
 640
 641 void ni_mio_ai_postgain_cal(void)
 642 {
 643         linear_fit_t l;
 644         double offset_r0;
 645         double offset_r7;
 646         double gain;
 647         double a;
 648
 649         check_gain_chan_x(&l,0,0,1);
 650         offset_r0=linear_fit_func_y(&l,caldacs[1].current);
 651         printf("offset r0 %g\n",offset_r0);
 652
 653         check_gain_chan_x(&l,0,7,1);
 654         offset_r7=linear_fit_func_y(&l,caldacs[1].current);
 655         printf("offset r7 %g\n",offset_r7);
 656
 657         gain=l.slope;
 658
 659         a=(offset_r0-offset_r7)/(200.0-1.0);
 660         a=caldacs[1].current-a/gain;
 661
 662         printf("%g\n",a);
 663
 664         caldacs[1].current=rint(a);
 665         update_caldac(1);
 666 }
 667
 668 void ni_mio_ai_postgain_cal_2(int chan,int dac,int range_lo,int range_hi,double gain)
 669 {
 670         double offset_lo,offset_hi;
 671         linear_fit_t l;
 672         double slope;
 673         double a;
 674
 675         check_gain_chan_x(&l,chan,range_lo,dac);
 676         offset_lo=linear_fit_func_y(&l,caldacs[dac].current);
 677         printf("offset lo %g\n",offset_lo);
 678
 679         check_gain_chan_x(&l,chan,range_hi,dac);
 680         offset_hi=linear_fit_func_y(&l,caldacs[dac].current);
 681         printf("offset hi %g\n",offset_hi);
 682
 683         slope=l.slope;
 684
 685         a=(offset_lo-offset_hi)/(gain-1.0);
 686         a=caldacs[dac].current-a/slope;
 687
 688         printf("%g\n",a);
 689
 690         caldacs[dac].current=rint(a);
 691         update_caldac(dac);
 692 }
 693
 694 void chan_cal(int adc,int cdac,int range,double target)
 695 {
 696         linear_fit_t l;
 697         double offset;
 698         double gain;
 699         double a;
 700
 701         check_gain_chan_x(&l,adc,range,cdac);
 702         offset=linear_fit_func_y(&l,caldacs[cdac].current);
 703         printf("offset %g\n",offset);
 704         gain=l.slope;
 705
 706         a=caldacs[cdac].current+(target-offset)/gain;
 707         printf("%g\n",a);
 708
 709         caldacs[cdac].current=rint(a);
 710         update_caldac(cdac);
 711 }
 712
 713
 714
 715 void channel_dependence(int adc,int range)
 716 {
 717         int i;
 718         double gain;
 719
 720         for(i=0;i<n_caldacs;i++){
 721                 gain=check_gain_chan(adc,range,i);
 722         }
 723 }
 724
 725 void caldac_dependence(int caldac)
 726 {
 727         int i;
 728         double gain;
 729
 730         for(i=0;i<8;i++){
 731                 gain=check_gain_chan(i,0,caldac);
 732         }
 733 }
 734
 735 void dump_curve(int adc,int caldac)
 736 {
 737         linear_fit_t l;
 738
 739         check_gain_chan_x(&l,adc,0,caldac);
 740 }
 741
 742
 743 void setup_caldacs(void)
 744 {
 745         int s,n_chan,i;
 746
 747         s=comedi_find_subdevice_by_type(dev,COMEDI_SUBD_CALIB,0);
 748         if(s<0){
 749                 printf("no calibration subdevice\n");
 750                 return;
 751         }
 752         //printf("calibration subdevice is %d\n",s);
 753
 754         n_chan=comedi_get_n_channels(dev,s);
 755
 756         for(i=0;i<n_chan;i++){
 757                 caldacs[i].subdev=s;
 758                 caldacs[i].chan=i;
 759                 caldacs[i].maxdata=comedi_get_maxdata(dev,s,i);
 760                 caldacs[i].current=0;
 761                 //printf("caldac %d, %d\n",i,caldacs[i].maxdata);
 762         }
 763
 764         n_caldacs=n_chan;
 765 }
 766
 767 void reset_caldacs(void)
 768 {
 769         int i;
 770
 771         for(i=0;i<n_caldacs;i++){
 772                 caldacs[i].current=caldacs[i].maxdata/2;
 773                 update_caldac(i);
 774         }
 775 }
 776
 777 void update_caldac(int i)
 778 {
 779         int ret;
 780
 781         //printf("update %d %d %d\n",caldacs[i].subdev,caldacs[i].chan,caldacs[i].current);
 782         //ret=comedi_data_write(dev,caldacs[i].subdev,caldacs[i].chan,0,0,caldacs[i].current);
 783         write_caldac(dev,caldacs[i].subdev,caldacs[i].chan,caldacs[i].current);
 784         ret=0;
 785         if(ret<0)
 786                 printf("comedi_data_write() error\n");
 787 }
 788
 789
 790 void check_gain(int ad_chan,int range)
 791 {
 792         int i;
 793
 794         for(i=0;i<n_caldacs;i++){
 795                 check_gain_chan(ad_chan,range,i);
 796         }
 797 }
 798
 799 double check_gain_chan(int ad_chan,int range,int cdac)
 800 {
 801         linear_fit_t l;
 802
 803         return check_gain_chan_x(&l,ad_chan,range,cdac);
 804 }
 805
 806 double check_gain_chan_x(linear_fit_t *l,int ad_chan,int range,int cdac)
 807 {
 808         int orig,i,n;
 809         int step;
 810         new_sv_t sv;
 811         double sum_err;
 812         int sum_err_count=0;
 813
 814         n=caldacs[cdac].maxdata+1;
 815         memset(l,0,sizeof(*l));
 816
 817         step=n/256;
 818         if(step<1)step=1;
 819         l->n=0;
 820
 821         l->y_data=malloc(n*sizeof(double)/step);
 822         if(l->y_data == NULL)
 823         {
 824                 perror("comedi_calibrate");
 825                 exit(1);
 826         }
 827
 828         orig=caldacs[cdac].current;
 829
 830         new_sv_init(&sv,dev,0,ad_chan,range,AREF_OTHER);
 831
 832         caldacs[cdac].current=0;
 833         update_caldac(cdac);
 834
 835         new_sv_measure(&sv);
 836         usleep(100000);
 837
 838         sum_err=0;
 839         for(i=0;i*step<n;i++){
 840                 caldacs[cdac].current=i*step;
 841                 update_caldac(cdac);
 842
 843                 new_sv_measure(&sv);
 844
 845                 l->y_data[i]=sv.average;
 846                 if(!isnan(sv.average)){
 847                         sum_err+=sv.error;
 848                         sum_err_count++;
 849                 }
 850                 l->n++;
 851         }
 852
 853         caldacs[cdac].current=orig;
 854         update_caldac(cdac);
 855
 856         l->yerr=sum_err/sum_err_count;
 857         l->dx=step;
 858         l->x0=0;
 859
 860         linear_fit_monotonic(l);
 861
 862         printf("caldac[%d] gain=%g V/bit err=%g S_min=%g dof=%g min=%g max=%g\n",
 863                 cdac,l->slope,l->err_slope,l->S_min,l->dof,l->min,l->max);
 864         //printf("--> %g\n",fabs(l.slope/l.err_slope));
 865
 866         if(verbose>=2){
 867                 static int dump_number=0;
 868                 double x,y;
 869
 870                 printf("start dump %d\n",dump_number);
 871                 for(i=0;i<l->n;i++){
 872                         x=l->x0+i*l->dx-l->ave_x;
 873                         y=l->y_data[i];
 874                         printf("D%d: %d %g %g %g\n",dump_number,i,y,
 875                                 l->ave_y+l->slope*x,
 876                                 l->ave_y+l->slope*x-y);
 877                 }
 878                 printf("end dump\n");
 879                 dump_number++;
 880         }
 881
 882
 883         free(l->y_data);
 884
 885         return l->slope;
 886 }
 887
 888
 889
 890
 891 /* helpers */
 892
 893
 894 int read_eeprom(int addr)
 895 {
 896         unsigned int data=0;
 897
 898         comedi_data_read(dev,eeprom_subdev,addr,0,0,&data);
 899
 900         return data;
 901 }
 902
 903 double read_chan(int adc,int range)
 904 {
 905         int n;
 906         new_sv_t sv;
 907
 908         new_sv_init(&sv,dev,0,adc,range,AREF_OTHER);
 909         sv.order=7;
 910         n=new_sv_measure(&sv);
 911
 912         printf("chan=%d ave=%g error=%g\n",adc,sv.average,sv.error);
 913
 914         return sv.average;
 915 }
 916
 917 void write_caldac(comedi_t *dev,int subdev,int addr,int val)
 918 {
 919         comedi_data_write(dev,subdev,addr,0,0,val);
 920 }
 921
 922
 923
 924 int new_sv_init(new_sv_t *sv,comedi_t *dev,int subdev,int chan,int range,int aref)
 925 {
 926         memset(sv,0,sizeof(*sv));
 927
 928         sv->subd=subdev;
 929         //sv->t.flags=TRIG_DITHER;
 930         sv->chan=chan;
 931         sv->range=range;
 932         sv->aref=aref;
 933
 934         //sv->chanlist[0]=CR_PACK(chan,range,aref);
 935
 936         sv->maxdata=comedi_get_maxdata(dev,subdev,chan);
 937         sv->rng=comedi_get_range(dev,subdev,chan,range);
 938
 939         sv->order=7;
 940
 941         return 0;
 942 }
 943
 944 int comedi_data_read_n(comedi_t *it,unsigned int subdev,unsigned int chan,
 945         unsigned int range,unsigned int aref,lsampl_t *data,unsigned int n)
 946 {
 947         comedi_insn insn;
 948         int ret;
 949
 950         if(n==0)return 0;
 951
 952         insn.insn = INSN_READ;
 953         insn.n = n;
 954         insn.data = data;
 955         insn.subdev = subdev;
 956         insn.chanspec = CR_PACK(chan,range,aref);
 957         /* enable dithering */
 958         insn.chanspec |= (1<<26);
 959
 960         ret = comedi_do_insn(it,&insn);
 961
 962         if(ret>0)return n;
 963
 964         printf("insn barfed: subdev=%d, chan=%d, range=%d, aref=%d, "
 965                 "n=%d, ret=%d, %s\n",subdev,chan,range,aref,n,ret,
 966                 strerror(errno));
 967
 968         return ret;
 969 }
 970
 971 int new_sv_measure(new_sv_t *sv)
 972 {
 973         lsampl_t *data;
 974         int n,i,ret;
 975
 976         double a,x,s,s2;
 977
 978         n=1<<sv->order;
 979
 980         data=malloc(sizeof(lsampl_t)*n);
 981         if(data == NULL)
 982         {
 983                 perror("comedi_calibrate");
 984                 exit(1);
 985         }
 986
 987         for(i=0;i<n;){
 988                 ret = comedi_data_read_n(dev,sv->subd,sv->chan,sv->range,
 989                         sv->aref,data+i,n-i);
 990                 if(ret<0){
 991                         printf("barf\n");
 992                         goto out;
 993                 }
 994                 i+=ret;
 995         }
 996
 997         s=0;
 998         s2=0;
 999         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
1000         for(i=0;i<n;i++){
1001                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
1002                 s+=x-a;
1003                 s2+=(x-a)*(x-a);
1004         }
1005         sv->average=a+s/n;
1006         sv->stddev=sqrt(n*s2-s*s)/n;
1007         sv->error=sv->stddev/sqrt(n);
1008
1009         ret=n;
1010
1011 out:
1012         free(data);
1013
1014         return ret;
1015 }
1016
1017 int new_sv_measure_order(new_sv_t *sv,int order)
1018 {
1019         lsampl_t *data;
1020         int n,i,ret;
1021         double a,x,s,s2;
1022
1023         n=1<<order;
1024
1025         data=malloc(sizeof(lsampl_t)*n);
1026         if(data == NULL)
1027         {
1028                 perror("comedi_calibrate");
1029                 exit(1);
1030         }
1031
1032         for(i=0;i<n;){
1033                 ret = comedi_data_read_n(dev,sv->subd,sv->chan,sv->range,
1034                         sv->aref,data+i,n-i);
1035                 if(ret<0){
1036                         printf("barf order\n");
1037                         goto out;
1038                 }
1039                 i+=ret;
1040         }
1041
1042         s=0;
1043         s2=0;
1044         a=comedi_to_phys(data[0],sv->rng,sv->maxdata);
1045         for(i=0;i<n;i++){
1046                 x=comedi_to_phys(data[i],sv->rng,sv->maxdata);
1047                 s+=x-a;
1048                 s2+=(x-a)*(x-a);
1049         }
1050         sv->average=a+s/n;
1051         sv->stddev=sqrt(n*s2-s*s)/n;
1052         sv->error=sv->stddev/sqrt(n);
1053
1054         ret=n;
1055
1056 out:
1057         free(data);
1058
1059         return ret;
1060 }
1061
1062
1063
1064 /* linear fitting */
1065
1066 int calculate_residuals(linear_fit_t *l);
1067
1068 int linear_fit_monotonic(linear_fit_t *l)
1069 {
1070         double x,y;
1071         double sxp;
1072         int i;
1073
1074         l->min=HUGE_VAL;
1075         l->max=-HUGE_VAL;
1076         l->s1=0;
1077         l->sx=0;
1078         l->sy=0;
1079         l->sxy=0;
1080         l->sxx=0;
1081         for(i=0;i<l->n;i++){
1082                 x=l->x0+i*l->dx;
1083                 y=l->y_data[i];
1084
1085                 if(isnan(y))continue;
1086
1087                 if(l->y_data[i]<l->min)l->min=l->y_data[i];
1088                 if(l->y_data[i]>l->max)l->max=l->y_data[i];
1089                 l->s1+=1;
1090                 l->sx+=x;
1091                 l->sy+=y;
1092                 l->sxy+=x*y;
1093                 l->sxx+=x*x;
1094         }
1095         sxp=l->sxx-l->sx*l->sx/l->s1;
1096
1097         l->ave_x=l->sx/l->s1;
1098         l->ave_y=l->sy/l->s1;
1099         l->slope=(l->s1*l->sxy-l->sx*l->sy)/(l->s1*l->sxx-l->sx*l->sx);
1100         l->err_slope=l->yerr/sqrt(sxp);
1101         l->err_ave_y=l->yerr/sqrt(l->s1);
1102
1103         calculate_residuals(l);
1104
1105         return 0;
1106 }
1107
1108 int calculate_residuals(linear_fit_t *l)
1109 {
1110         double x,y;
1111         double res,sum_res2;
1112         int i;
1113
1114         sum_res2=0;
1115         for(i=0;i<l->n;i++){
1116                 x=l->x0+i*l->dx-l->ave_x;
1117                 y=l->y_data[i];
1118
1119                 if(isnan(y))continue;
1120
1121                 res=l->ave_y+l->slope*x-y;
1122
1123                 sum_res2+=res*res;
1124         }
1125         l->S_min=sum_res2/(l->yerr*l->yerr);
1126         l->dof=l->s1-2;
1127
1128         return 0;
1129 }
1130
1131 double linear_fit_func_y(linear_fit_t *l,double x)
1132 {
1133         return l->ave_y+l->slope*(x-l->ave_x);
1134 }
1135