ref: 70c85915eace83142b84e4f65f5db421cf0c09e3
dir: /test/corr.c/
/*
corr.c
print correlation coeff's for 2 input files of 'short' values
Copyright (C) 1999 Stanley J. Brooks <stabro@megsinet.net>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <string.h>
// for open,read,write:
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <errno.h>
#include <unistd.h>
// for clock/time:
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define BLEN 0x20000
#define HLEN 0x1000
static u_int32_t dist[0x10000];
#define dis0 (dist+0x8000)
static short buff1[BLEN];
static short buff2[BLEN];
static short hist[HLEN];
inline static double square(double x) {return (x*x);}
static void readin(int fd, short *p, int bs)
{
int r;
do {
r=read(fd,(char*)p,2*bs);
}while (r==-1 && errno==EINTR);
if (r==-1 || r!=2*bs) {
perror("Input error");
exit(1);
}
}
int main(int argct, char **argv)
{
int fd1,fd2;
int avgn=0,assay=0,code=0;
int hx=0,hsum=0;
int len,len1,len2,x,r;
char *fnam1,*fnam2;
double v1=0,v2=0,s1=0;
/*
* Parse the command line arguments
*/
while (*++argv && **argv == '-')
while (*(++*argv))
switch (**argv) {
case 'a':
assay=1;
break;
case 'c':
code=1;
break;
case 'h':
fprintf(stderr,"./stat [-n avgn] <file>\n");
break;
case 'n':
if (!*++*argv) { // end of this param string
if (!argv[1]) break;
++argv;
}
avgn=atoi(*argv);
*argv += strlen(*argv)-1; // skip to end-1 of this param string
break;
}
if (avgn<=0) avgn=1; else if (avgn>HLEN) avgn=HLEN;
bzero(hist,sizeof(hist));
fnam1=*argv++;
fd1=open(fnam1,O_RDWR);
if (fd1<0) {
fprintf(stderr,"Open: %s %s\n",fnam1,strerror(errno)); return(1);
}
len1=lseek(fd1,0,SEEK_END);
r=lseek(fd1,0,SEEK_SET);
fnam2=*argv++;
fd2=open(fnam2,O_RDWR);
if (fd2<0) {
fprintf(stderr,"Open: %s %s\n",fnam2,strerror(errno)); return(1);
}
len2=lseek(fd2,0,SEEK_END);
r=lseek(fd2,0,SEEK_SET);
bzero(dist,sizeof(dist));
len = (len1<len2)? len1:len2;
len /= 2; /* now len is number of shorts data to read */
for (x=0; x<len; ){
int bs;
int j;
double c11=0, c12=0, c22=0;
bs=len-x; if (bs>BLEN) bs=BLEN; /* number of shorts to read */
readin(fd1,buff1,bs);
readin(fd2,buff2,bs);
for (j=0; j<bs; j++) {
c11 += buff1[j]*buff1[j];
c12 += buff1[j]*buff2[j];
c22 += buff2[j]*buff2[j];
}
c11 /= bs;
c12 /= bs;
c22 /= bs;
{
double d11=(c11+2*c12+c22)/4;
double d22=(c11-2*c12+c22)/4;
double d12=(c11-c22)/4;
printf("%8.1f%8.1f cf=%f",sqrt(c11),sqrt(c22),c12/sqrt(c11*c22));
printf(" | %8.1f%8.1f cf=%f\n",sqrt(d11),sqrt(d22),d12/sqrt(d11*d22));
}
for (j=0; j<bs; j++) {
int y;
//y=(abs(buff1[j])<abs(buff2[j]))? buff1[j]:buff2[j];
y=(buff1[j]+buff2[j]+1)/2;
hsum -= hist[hx];
hsum += y;
hist[hx] = y;
if (++hx == avgn) hx=0;
s1 += abs(y);
v1 += hsum;
v2 += square(hsum);
y = (hsum+avgn/2)/avgn;
dis0[y]++;
}
x += bs;
}
v1 /= len;
v2 /= len;
printf("%8d %5d=n %10.2f=avg %10.2f=rms\n",len,avgn,v1/avgn,sqrt(v2-v1*v1)/avgn);
if (code) {
int32_t j,tot,cut;
cut = len/2;
for (tot=0,j=0; j<0x8000; j++) {
int32_t tot1;
tot1 = tot+dis0[j];
if (j!=0) tot1 += dis0[-j];
if (tot1>cut && tot<=cut) {
printf(" |%d| %d of %d\n",j-1,tot,cut);
cut += (len-cut)/2;
}
tot=tot1;
}
}
if (assay)
for (x=-0x8000; x<=0x7fff; x++) {
printf("%6d %6d\n",x,dis0[x]);
}
return 0;
}