Low level routine
fr -
eng
pspect_c - inter/cross correlation spectra computational routine
/* pspect_c subroutine
*
* 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"
/* pspect_c estimation de la densité spectrale de puissance par la méthode
* du périodogramme.
*
* sec_l : taille de la fenetre
* n_sec : nombre de moyennes désirées
* r_io : adresse de départ de la partie réelle du vecteur temporel d'entrée, longueur sec_l
* i_io : adresse de départ de la partie imaginaire du vecteur temporel d'entrée, longueur sec_l
* r_out : tableau de sortie (état), longueur sec_l
* n_call : numéro d'appel courant (état)
* win : tableau des coef de la fenetre de ponderation, longueur sec_l
* fft_pr : structure pour calcul de fft
*
* rmq : seul la partie réelle est considérée ici.
*/
void pspect_c(int sec_l,int n_sec,fft_pr_struct *fft_pr, \
int *n_call,double *win,double *r_io,double *i_io,double *r_out)
{
/* define local variable */
double sum;
int i;
/* sum entries */
sum = 0.;
for(i=0;i<sec_l;i++) {
sum += r_io[i];
}
sum = sum/sec_l;
/* normalize and mult by coef of win */
for(i=0;i<sec_l;i++) {
r_io[i] = win[i]*(r_io[i]-sum);
}
/* execute fft computation */
/* calc_fft(sec_l,isn,ffttyp,n_work,work->work,r_io, i_io);*/
calc_fft_c(fft_pr,r_io,i_io);
/* compute conjugate of the complex spectrum
* and add it to the previous spectrum
*/
for(i=0;i<sec_l;i++) {
r_out[i] = r_out[i]+(r_io[i]*r_io[i]+i_io[i]*i_io[i]);
}
/* update average counter */
(*n_call)++;
/* if average counter cross user average param
* then display spectrum
*/
if ((*n_call) == n_sec) {
/* compute power of window */
sum = 0.;
for(i=0;i<sec_l;i++) {
sum += win[i];
}
sum = sum * (*n_call);
/* normalize */
for(i=0;i<sec_l;i++) {
r_out[i] = r_out[i]/sum;
}
/* reinit average counter */
(*n_call) = 0;
}
}
A. Layec