Routine bas-niveau
fr - eng

pspect_c - routine de calcul de spectre par inter/cross corélation

Module

• Modnum

Contenu du fichier

/* 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;
  }
}

Auteurs

A. Layec A. Layec