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
square(double x)46 inline static double square(double x) {return (x*x);}
47
readin(int fd,short * p,int bs)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
main(int argct,char ** argv)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
188