test/corr.c

   1 /*
   2   corr.c
   3     print correlation coeff's for 2 input files of 'short' values
   4
   5         Copyright (C) 1999 Stanley J. Brooks <stabro@megsinet.net>
   6
   7   This program is free software; you can redistribute it and/or modify
   8   it under the terms of the GNU General Public License as published by
   9   the Free Software Foundation; either version 2 of the License, or
  10   (at your option) any later version.
  11
  12   This program is distributed in the hope that it will be useful,
  13   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15   GNU General Public License for more details.
  16
  17   You should have received a copy of the GNU General Public License
  18   along with this program; if not, write to the Free Software
  19   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20
  21 */
  22
  23 #include <string.h>
  24 // for open,read,write:
  25 #include <sys/types.h>
  26 #include <sys/stat.h>
  27 #include <fcntl.h>
  28 #include <errno.h>
  29 #include <unistd.h>
  30 // for clock/time:
  31 #include <time.h>
  32 #include <stdlib.h>
  33 #include <stdio.h>
  34 #include <math.h>
  35
  36 #define BLEN 0x20000
  37 #define HLEN 0x1000
  38
  39 static u_int32_t dist[0x10000];
  40 #define dis0 (dist+0x8000)
  41
  42 static short buff1[BLEN];
  43 static short buff2[BLEN];
  44 static short hist[HLEN];
  45
  46 inline static double square(double x) {return (x*x);}
  47
  48 static void readin(int fd, short *p, int bs)
  49 {
  50         int r;
  51         do {
  52                 r=read(fd,(char*)p,2*bs);
  53         }while (r==-1 && errno==EINTR);
  54         if (r==-1 || r!=2*bs) {
  55                 perror("Input error");
  56                 exit(1);
  57         }
  58 }
  59
  60 int main(int argct, char **argv)
  61 {
  62         int fd1,fd2;
  63         int avgn=0,assay=0,code=0;
  64         int hx=0,hsum=0;
  65         int len,len1,len2,x,r;
  66         char *fnam1,*fnam2;
  67         double v1=0,v2=0,s1=0;
  68
  69         /*
  70          * Parse the command line arguments
  71          */
  72         while (*++argv && **argv == '-')
  73                 while (*(++*argv))
  74                         switch (**argv) {
  75                         case 'a':
  76                                 assay=1;
  77                                 break;
  78                         case 'c':
  79                                 code=1;
  80                                 break;
  81                         case 'h':
  82                                 fprintf(stderr,"./stat [-n avgn] <file>\n");
  83                                 break;
  84                         case 'n':
  85                                 if (!*++*argv) { // end of this param string
  86                                         if (!argv[1]) break;
  87                                         ++argv;
  88                                 }
  89                                 avgn=atoi(*argv);
  90                                 *argv += strlen(*argv)-1;       // skip to end-1 of this param string
  91                                 break;
  92                         }
  93
  94         if (avgn<=0) avgn=1; else if (avgn>HLEN) avgn=HLEN;
  95         bzero(hist,sizeof(hist));
  96
  97         fnam1=*argv++;
  98         fd1=open(fnam1,O_RDWR);
  99         if (fd1<0) {
 100                 fprintf(stderr,"Open: %s %s\n",fnam1,strerror(errno)); return(1);
 101         }
 102         len1=lseek(fd1,0,SEEK_END);
 103         r=lseek(fd1,0,SEEK_SET);
 104
 105         fnam2=*argv++;
 106         fd2=open(fnam2,O_RDWR);
 107         if (fd2<0) {
 108                 fprintf(stderr,"Open: %s %s\n",fnam2,strerror(errno)); return(1);
 109         }
 110         len2=lseek(fd2,0,SEEK_END);
 111         r=lseek(fd2,0,SEEK_SET);
 112
 113         bzero(dist,sizeof(dist));
 114         len = (len1<len2)? len1:len2;
 115         len /= 2; /* now len is number of shorts data to read */
 116         for (x=0; x<len; ){
 117                 int bs;
 118                 int j;
 119                 double c11=0, c12=0, c22=0;
 120
 121                 bs=len-x; if (bs>BLEN) bs=BLEN; /* number of shorts to read */
 122
 123                 readin(fd1,buff1,bs);
 124                 readin(fd2,buff2,bs);
 125
 126                 for (j=0; j<bs; j++) {
 127                         c11 += buff1[j]*buff1[j];
 128                         c12 += buff1[j]*buff2[j];
 129                         c22 += buff2[j]*buff2[j];
 130                 }
 131                 c11 /= bs;
 132                 c12 /= bs;
 133                 c22 /= bs;
 134
 135                 {
 136                         double d11=(c11+2*c12+c22)/4;
 137                         double d22=(c11-2*c12+c22)/4;
 138                         double d12=(c11-c22)/4;
 139
 140                         printf("%8.1f%8.1f cf=%f",sqrt(c11),sqrt(c22),c12/sqrt(c11*c22));
 141                         printf(" | %8.1f%8.1f cf=%f\n",sqrt(d11),sqrt(d22),d12/sqrt(d11*d22));
 142                 }
 143
 144                 for (j=0; j<bs; j++) {
 145                         int y;
 146                         //y=(abs(buff1[j])<abs(buff2[j]))? buff1[j]:buff2[j];
 147                         y=(buff1[j]+buff2[j]+1)/2;
 148                         hsum -= hist[hx];
 149                         hsum += y;
 150                         hist[hx] = y;
 151                         if (++hx == avgn) hx=0;
 152                         s1 += abs(y);
 153                         v1 += hsum;
 154                         v2 += square(hsum);
 155                         y = (hsum+avgn/2)/avgn;
 156                         dis0[y]++;
 157                 }
 158
 159                 x += bs;
 160         }
 161         v1 /= len;
 162         v2 /= len;
 163         printf("%8d %5d=n %10.2f=avg %10.2f=rms\n",len,avgn,v1/avgn,sqrt(v2-v1*v1)/avgn);
 164
 165         if (code) {
 166                 int32_t j,tot,cut;
 167                 cut = len/2;
 168                 for (tot=0,j=0; j<0x8000; j++) {
 169                         int32_t tot1;
 170                         tot1 = tot+dis0[j];
 171                         if (j!=0) tot1 += dis0[-j];
 172                         if (tot1>cut && tot<=cut) {
 173                                 printf(" |%d| %d of %d\n",j-1,tot,cut);
 174                                 cut += (len-cut)/2;
 175                         }
 176                         tot=tot1;
 177                 }
 178         }
 179
 180         if (assay)
 181                 for (x=-0x8000; x<=0x7fff; x++) {
 182                         printf("%6d %6d\n",x,dis0[x]);
 183                 }
 184
 185         return 0;
 186 }
 187