/* corr_c subroutine * discrete cross-correlation computation * function * * Copyright (C) 2007-2009 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; }