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