Low level routine
fr - eng


corr_c - cross-correlation computational routine

Module

File content


/* corr_c subroutine
 * discrete cross-correlation computation
 * function
 *
 * Copyright (C) 2007-2011 Alan Layec
 *
 * This file is part of modnumlib.
 *
 * modnumlib 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.
 *
 * modnumlib 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 modnumlib; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 */

/* REVISION HISTORY :
 * $Log$
 */

#include "modnum_lib.h"

/*
 * corr_c routine de calcul d'une intercorrelation discrète
 *
 * Entrées :
 * n       : dimension 1 des matrices en entrée (entier,scalaire)
 * m       : dimension 2 des matrices en entrée (entier,scalaire)
 * typ     : 0 : les moyennes sont passées en argument
 *           1 : les moyennes doivent être calculées
 * lag     : longueur de l'intercorrelation (entier,scalaire)
 * mean_u1 : moyenne de la série u1 (vecteur double de taille m) (E/S)
 * mean_u2 : moyenne de la série u2 (vetceur double de taille m) (E/S)
 * u1      : matrice 1 en entrée (double)
 * u2      : matrice 2 en entrée (double)
 *
 * Sorties :
 * y   : matrice en sortie (double) de taille lag,m
 *
 */
void corr_c(int *n,int *m,int *typ,int *lag,double *mean_u1,double *mean_u2,\
            double *u1,double *u2,double *y)
{
 /*Déclaration des variables compteurs*/
 int i,j,l;

 int ind_u;
 int ind_y;

 /*Déclaration de variables auxiliares*/
 double y_tmp; /* valeur de l'intercorr sur un indice donné */
 double N_1 = 1/((double)(*n)); /* facteur d'equiprobabilité */

 switch (*typ)
 {
   /* cas moyennes données */
   case 0 : {
             /* calcul intercorr */
             for(l=0;l<(*m);l++) {

               ind_u = (*n)*l;
               ind_y = (*lag)*l;

               for(i=0;i<(*lag);i++) {
                 /* initialise y_tmp */
                 y_tmp = 0.;

                 /* calcul intercorr de l'indice i */
                 for(j=0;j<((*n)-i);j++) {
                   y_tmp += (u1[ind_u+j] - mean_u1[l])*(u2[ind_u+j+i] - mean_u2[l]);
                 }

                 /* calcul la sortie normalisée */
                 y[ind_y+i] = N_1 * y_tmp;
               }
             }
             break;
            }

   /* cas moyennes à calculer */
   case 1 : {
             for(l=0;l<(*m);l++) {

               ind_u = (*n)*l;
               ind_y = (*lag)*l;

               /* init moyennes des séries m */
               mean_u1[l] = 0.;
               mean_u2[l] = 0.;

               /* calcul moyennes */
               for(i=0;i<(*n);i++) {
                 mean_u1[l] += u1[ind_u+i];
                 mean_u2[l] += u2[ind_u+i];
               }

               /* normalisation des moyennes */
               mean_u1[l] *= N_1;
               mean_u2[l] *= N_1;


               /* calcul intercorr */
               for(i=0;i<(*lag);i++) {

                 /* initialise y_tmp */
                 y_tmp = 0.;

                 /* calcul intercorr de l'indice i */
                 for(j=0;j<((*n)-i);j++) {
                   y_tmp += (u1[ind_u+j] - mean_u1[l])*(u2[ind_u+j+i] - mean_u2[l]);
                 }

                 /* calcul la sortie normalisée */
                 y[ind_y+i] = y_tmp * N_1;
               }
             }
             break;
            }

 }

 return;
}

void corri_c(int *n,int *m,int *typ,int *lag,double *mean_u1,double *mean_u2,\
             int *u1,int *u2,double *y)
{
 /*Déclaration des variables compteurs*/
 int i,j,l;

 int ind_u;
 int ind_y;

 /*Déclaration de variables auxiliares*/
 double y_tmp; /* valeur de l'intercorr sur un indice donné */
 double N_1 = 1/((double)(*n)); /* facteur d'equiprobabilité */

 switch (*typ)
 {
   /* cas moyennes données */
   case 0 : {
             /* calcul intercorr */
             for(l=0;l<(*m);l++) {

               ind_u = (*n)*l;
               ind_y = (*lag)*l;

               for(i=0;i<(*lag);i++) {
                 /* initialise y_tmp */
                 y_tmp = 0.;

                 /* calcul intercorr de l'indice i */
                 for(j=0;j<((*n)-i);j++) {
                   y_tmp += ((double)u1[ind_u+j] - mean_u1[l])*\
                            ((double)u2[ind_u+j+i] - mean_u2[l]);
                 }

                 /* calcul la sortie normalisée */
                 y[ind_y+i] = N_1 * y_tmp;
               }
             }
             break;
            }

   /* cas moyennes à calculer */
   case 1 : {
             for(l=0;l<(*m);l++) {

               ind_u = (*n)*l;
               ind_y = (*lag)*l;

               /* init moyennes des séries m */
               mean_u1[l] = 0.;
               mean_u2[l] = 0.;

               /* calcul moyennes */
               for(i=0;i<(*n);i++) {
                 mean_u1[l] += (double)u1[ind_u+i];
                 mean_u2[l] += (double)u2[ind_u+i];
               }

               /* normalisation des moyennes */
               mean_u1[l] *= N_1;
               mean_u2[l] *= N_1;

               /* calcul intercorr */
               for(i=0;i<(*lag);i++) {
                 /* initialise y_tmp */
                 y_tmp = 0.;

                 /* calcul intercorr de l'indice i */
                 for(j=0;j<((*n)-i);j++) {
                   y_tmp += ((double)u1[ind_u+j] - mean_u1[l])*\
                                    ((double)u2[ind_u+j+i] - mean_u2[l]);
                 }

                 /* calcul la sortie normalisée */
                 y[ind_y+i] = y_tmp * N_1;
               }
             }
             break;
            }

 }

 return;
}

Authors

A. Layec