home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_09_05 / 9n05059a < prev    next >
Text File  |  1991-02-07  |  11KB  |  508 lines
  1. /************************************* 
  2.  *                                   * 
  3.  *    LINEAR & NONLINEAR VARIANCE    * 
  4.  *            TEST PROGRAM           * 
  5.  *                                   * 
  6.  *  M.E. Brandt, Ph.D.               * 
  7.  *  02/06/91 Revision                * 
  8.  *                                   * 
  9.  *************************************/
  10.  
  11. #include <stdio.h>
  12. #include <math.h>
  13. #include <malloc.h>
  14.  
  15.  
  16. #define NMAX                   1024
  17. #define NBINS                    10
  18. #define MAX_DELAY                30
  19. #define MAXRAND             32767.0
  20. #define DBL_LEN      sizeof(double)
  21. #define PI        3.141592653589793
  22. #define TWOPI            (2.0 * PI)
  23.  
  24.  
  25.  
  26. double linvar(double *y, double *x, int nn);
  27. double nonlinvar(double *x, double *y, int n, int nbins);
  28. int    number_bins(int x);
  29. void   estimate_delay(double *x, double *y, int n);
  30.  
  31.  
  32.  
  33. main()
  34. {
  35.    
  36.    int xx, i, j,  nb, n, ch;
  37.  
  38.    double *x, *y, scale, r2, h2, start, 
  39.           noise, gain, a, inc, mean, num;
  40.       
  41.    char u[3], c;
  42.    
  43.    
  44.  
  45.    srand(1);
  46.    
  47.         
  48.   /* SELECT A FUNCTION TO TEST */
  49.   
  50.   do {
  51.   
  52.         printf("\n\nSelect a y vs. x function and ");
  53.         printf("compute variances:\n");
  54.         printf("\n1. y = ax + noise (linear)");
  55.         printf("\n2. y = sin(x) + noise");
  56.         printf("\n3. y = sgn(x)*(x*x) + noise");
  57.         printf("\n4. y = pow(3*x, 3) + noise");
  58.         printf("\n5. y = pow(5*x, 5) + noise");
  59.         printf("\n6. Exit");
  60.         printf("\n\nSelect 1-6: ");
  61.         scanf("%d", &ch);
  62.         
  63.         if(ch == 6) exit(0);
  64.  
  65.         do {
  66.                 
  67.             printf("\n\nHow many data points?: ");
  68.             scanf("%d", &n);
  69.         }
  70.         while(n > NMAX);
  71.         
  72.         do {
  73.         
  74.             printf("\n\nEnter noise gain: ");
  75.             scanf("%lf", &gain);
  76.          }
  77.          while(gain < 0.0 || gain > 1.0);
  78.                 
  79.         /* allocate space for x and y vectors */
  80.         
  81.         if((x = (double *)malloc(n * DBL_LEN)) == NULL) {
  82.             fprintf(stderr, "\nMalloc error; x\n");
  83.             exit(1);
  84.         }
  85.         
  86.         if((y = (double *)malloc(n * DBL_LEN)) == NULL) {
  87.             fprintf(stderr, "\nMalloc error; y\n");
  88.             exit(1);
  89.         }
  90.  
  91.  
  92.       /* generate a random uniformly distributed x[n] 
  93.          between + or - 0.5 
  94.       */
  95.     
  96.       for(i=0; i<n; i++) x[i] = 
  97.                           ((double)rand()/MAXRAND - 0.5);
  98.     
  99.                   
  100.       /* choose y[n] */
  101.       
  102.       switch(ch) {
  103.       
  104.           case 1:  printf("\nEnter a: ");
  105.                    scanf("%lf", &a);
  106.                     
  107.                    for(i=0; i<n; i++)  {
  108.                        noise = gain * 
  109.                            ((double)rand()/MAXRAND - 0.5);
  110.                        y[i] = a * x[i] + noise;       
  111.                    }
  112.                    break;
  113.       
  114.           case 2:  for(i=0; i<n; i++)  {
  115.                        noise = gain * 
  116.                            ((double)rand()/MAXRAND - 0.5);
  117.                        y[i] = sin(TWOPI * x[i]) + noise;       
  118.                    }
  119.                    break;
  120.  
  121.           case 3:  for(i=0; i<n; i++)  {
  122.                        noise = gain * 
  123.                            ((double)rand()/MAXRAND - 0.5);
  124.                        y[i] = (x[i] * x[i]);
  125.                        if(x[i] < 0.0) y[i] *= -1.0;
  126.                        y[i] += noise;
  127.                     }
  128.                     break;
  129.     
  130.            case 4:  for(i=0; i<n; i++)  {
  131.                        noise = gain * 
  132.                            ((double)rand()/MAXRAND - 0.5);
  133.                        y[i] = pow(3.0*x[i], 3.0) + noise;
  134.                     }
  135.                     break;
  136.    
  137.            case 5:  for(i=0; i<n; i++)  {
  138.                        noise = gain * 
  139.                            ((double)rand()/MAXRAND - 0.5);
  140.                        y[i] = pow(5.0*x[i], 5.0) + noise;
  141.                     }
  142.                     break;
  143.      
  144.           default:  break;
  145.     
  146.         }
  147.  
  148.     printf("\nTest delay estimation <y/n>?: ");
  149.     scanf("%s", u);
  150.     c = toupper(*u);
  151.     if(c == 'Y') {
  152.        estimate_delay(x, y, n);
  153.        continue;
  154.     }
  155.  
  156.  
  157.             
  158.     /* compute linear variance of x vs. y */
  159.             
  160.     r2 = linvar(x, y, n);
  161.     
  162.     printf("\nLinear r2(x,y) = %f\n", r2);
  163.     
  164.     
  165.     /* choose number of bins */
  166.     
  167.     nb = number_bins(n);       
  168.       
  169.     /* compute nonlinear variance of x vs. y */
  170.     
  171.     h2 = nonlinvar(x, y, n, nb);
  172.         
  173.     printf("\nNonlinear h2(x,y) = %f\n", h2);
  174.  
  175.  
  176.     /* y vs. x */
  177.     
  178.     r2 = linvar(y, x, n);
  179.     
  180.     printf("\nLinear r2(y,x) = %f\n", r2);
  181.     
  182.     
  183.     h2 = nonlinvar(y, x, n, nb);
  184.         
  185.     printf("\nNonlinear h2(y,x) = %f\n", h2);
  186.  
  187.     free(x); free(y);  
  188.        
  189.   }
  190.   while(1);
  191.  
  192.   
  193. }
  194.         
  195.  
  196.  
  197.  
  198.  
  199.  
  200.  
  201. /* routine to find variances & delay between 2 lagged 
  202.    signals */
  203.    
  204. double buf[NMAX+MAX_DELAY], r2[MAX_DELAY+10],
  205.        h2[MAX_DELAY+10];
  206.        
  207. void estimate_delay(double *x, double *y, int n)
  208. {
  209.  
  210.    int i, j, d, delay, lag, max_lag, np, nb;   
  211.    
  212.    double max;
  213.  
  214.  
  215.    do {
  216.         
  217.       printf("\nEnter delay between x and y in samples: ");
  218.       scanf("%d", &delay);
  219.    }
  220.    while((d = abs(delay)) > MAX_DELAY);
  221.    
  222.    
  223.    for(i=0; i<d; i++) 
  224.        buf[i] = ((double)rand()/MAXRAND - 0.5);
  225.  
  226.    max_lag = d + 10;
  227.       
  228.    if(delay > 0) {   /* y lags x */
  229.  
  230.       /* move y to buf */
  231.          
  232.       for(j=0, i=d; j<n; i++, j++) buf[i] = y[j];
  233.  
  234.       for(lag=0; lag<max_lag; lag++) {
  235.    
  236.           np = n - lag;     /* compute across less points */
  237.         
  238.           /* get linear variance at lag */
  239.           
  240.           r2[lag] = linvar(x, &buf[lag], np);
  241.              
  242.           nb = number_bins(np);
  243.         
  244.           /* get nonlinear variance at lag */
  245.             
  246.           h2[lag] = nonlinvar(x, &buf[lag], np, nb);
  247.       
  248.       }
  249.       
  250.    }
  251.                    
  252.    else if(delay < 0) {  /* x lags y */
  253.    
  254.            for(j=0, i=d; j<n; i++, j++) buf[i] = x[j];
  255.  
  256.            for(lag=0; lag<max_lag; lag++) {
  257.    
  258.                np = n - lag;
  259.         
  260.                r2[lag] = linvar(&buf[lag], y, np);
  261.              
  262.                nb = number_bins(np);
  263.         
  264.                h2[lag] = nonlinvar(&buf[lag], y, np, nb);
  265.       
  266.            }
  267.  
  268.     }
  269.      
  270.  
  271.    /* find maximum r2 and corresponding delay */
  272.    
  273.       max = r2[0];
  274.       for(i=1; i<max_lag; i++) 
  275.           if(r2[i] > max) {max = r2[i]; j = i;} 
  276.  
  277.       if(delay < 0) j = -j;
  278.          
  279.       printf("\nTrue delay = %d; maximum r2 = %f \
  280. found @ sample %d\n",
  281.              delay, max, j);
  282.    
  283.    
  284.    /* find maximum h2 and corresponding delay */
  285.       
  286.       max = h2[0];
  287.       for(i=1; i<max_lag; i++) 
  288.           if(h2[i] > max) {max = h2[i]; j = i;} 
  289.    
  290.       if(delay < 0) j = -j;
  291.       
  292.       printf("\nTrue delay = %d; maximum h2 = %f \
  293. found @ sample %d\n",
  294.              delay, max, j);
  295.        
  296.  
  297. }
  298.  
  299.  
  300.  
  301. /* compute number of bins for nonlinear variance calc. */
  302.  
  303. int number_bins(int x)
  304. {
  305.    
  306.     int y;
  307.       
  308.     if(x < 129) y = NBINS - 4;
  309.     else
  310.     if(x < 257) y = NBINS - 3;
  311.     else
  312.     if(x < 513) y = NBINS - 2;
  313.     else        y = NBINS;
  314.      
  315.     return y;
  316.  
  317. }
  318.  
  319.     
  320.     
  321. /* routine to compute linear variance measure */
  322.  
  323. double linvar(double *y, double *x, int nn)
  324. {
  325.  
  326.    register int j;
  327.  
  328.    double z, smx, smy, sqx, sqy, sxy, corr;
  329.    double sqrt(), varx, vary, cov, nnn;
  330.  
  331.  
  332.  
  333. /* initialize values */
  334.  
  335. nnn = (double)nn;
  336.  
  337. corr = sqy = 0.0;
  338.  
  339.  
  340. smx = 0.0;
  341. smy = 0.0;
  342. sqx = 0.0;
  343. sxy = 0.0;
  344.  
  345. for(j=0; j<nn; j++)  {
  346.  
  347.    smx += x[j];
  348.    smy += y[j];
  349.    sqx += (x[j] * x[j]);
  350.    sxy += (x[j] * y[j]);
  351.  
  352.  
  353.    sqy += (y[j] * y[j]);
  354.  
  355.  
  356. } /* end for */
  357.  
  358. /* Compute covariance and variances */
  359.  
  360. cov  = (nnn * sxy) - (smx * smy);
  361. varx = (nnn * sqx) - (smx * smx);
  362.  
  363.  
  364.  
  365. vary = (nnn * sqy) - (smy * smy);
  366.  
  367. z = sqrt(varx * vary);
  368.  
  369. if(z != 0.0) corr = cov / z;
  370. else corr = 0.0;
  371.  
  372. return(corr*corr);
  373.  
  374. }
  375.  
  376.  
  377.  
  378. /* routine to compute nonlinear variance measure */
  379.  
  380. int bin[NBINS][2*NMAX/NBINS]; 
  381. double q[NBINS], p[NBINS];
  382.  
  383. double nonlinvar(double *x, double *y, int n, int nbins)
  384. {
  385.  
  386.     int i, j, k, l, index, index2, last;
  387.      
  388.     double min, max, range, width, hwidth, start, end,   
  389.            mean, mean2, totvar, unvar, s, f, h2, uu, offset;
  390.            
  391.            
  392.      /* find max and min of x array */
  393.      
  394.      min = max = x[0]; last = 0;
  395.      
  396.      for(i=1; i<n; i++) {
  397.          if(x[i] > max) {max = x[i]; last = i;}
  398.          else
  399.          if(x[i] < min) min = x[i];
  400.      }
  401.      
  402.      range = max - min;
  403.           
  404.      width = range/(double)nbins; 
  405.      
  406.      hwidth = width/2.0;
  407.      
  408.      for(i=0; i<nbins; i++) bin[i][0] = 0;
  409.              
  410.     /* get histogram of indexes */
  411.  
  412.     for(i=0; i<n; i++) {
  413.         for(j=0; j<nbins; j++) {
  414.             start = (double)j*width + min;
  415.             end = start + width;
  416.             
  417.             if((x[i] >= start) && (x[i] < end)) {
  418.                 ++bin[j][0]; 
  419.                 bin[j][bin[j][0]] = i;
  420.                 break;
  421.             }
  422.          }
  423.       }
  424.       
  425.       /* maximum x value (last one) */
  426.       
  427.       j = nbins-1;
  428.       ++bin[j][0];  
  429.       bin[j][bin[j][0]] = last;
  430.       
  431.       
  432.     /* compute x-midpoints and y average amplitudes */
  433.  
  434.      mean2 = 0.0;     
  435.      offset = hwidth + min;
  436.      
  437.      for(i=0; i<nbins; i++) {
  438.      
  439.          /* x value midpoints for each bin */
  440.              
  441.          p[i] =  ((double)i * width) + offset;
  442.          
  443.          /* corresponding y average amplitudes */
  444.          
  445.          mean = 0.0;
  446.          
  447.          for(k=1; k<=bin[i][0]; k++) mean += y[bin[i][k]];
  448.          
  449.          q[i] =  mean/(double)bin[i][0];
  450.          mean2 += mean;
  451.                  
  452.       }
  453.       
  454.             
  455.       mean = mean2/(double)n;    /* overall */
  456.       
  457.       
  458.             
  459.       /* compute h*h coefficient */
  460.       
  461.       /* first compute total variance of y */
  462.       
  463.       totvar = 0.0;
  464.       for(i=0; i<n; i++) {
  465.           s = y[i] - mean;
  466.           totvar += (s * s);
  467.       }
  468.  
  469.       /* compute total unexplained variance of y */
  470.       
  471.       unvar = 0.0;      
  472.       for(i=0; i<n; i++) {
  473.       
  474.           /* find f(x[i]), the nonlinear value */
  475.           
  476.           for(j=1; j<nbins-1; j++) {
  477.               if(x[i] <= p[j]) {
  478.                  index = j-1;
  479.                  goto out1;
  480.               }
  481.           }
  482.           
  483.           index = nbins-2;
  484.           
  485.    out1:       
  486.           index2 = index++;
  487.           
  488.           f = ((q[index2] - q[index])/
  489.                (p[index2] - p[index])) *
  490.                  (x[i] - p[index]) + q[index];
  491.           
  492.           uu = y[i] - f;
  493.           unvar += (uu * uu);         
  494.           
  495.      }
  496.                               
  497.      /* now compute nonlinear regression coefficient */
  498.      
  499.           
  500.      h2 = 1.0 - (unvar/totvar);
  501.      
  502.      return h2;
  503.           
  504. }
  505.            
  506.  
  507.              
  508.