source: Sophya/trunk/SigPredictor/Tabulatedfilter.cc

                        // Dominique YVON, CEA/DAPNIA/SPP 02/2000

#include "macenvvariables.h"
#include "tabulatedfilter.h"

NumRecipes  TabulatedFilter::NR;
static unsigned long LastfreqIndex;

TabulatedFilter::TabulatedFilter(char filename[])
{       LastfreqIndex=0;
        FILE* pfile;
//      char fileCur[150];
        
        // On lit le tableau des abscisses
        pfile=fopen(filename,"r");
        if (pfile==NULL)
        {       cerr<< "Erreur a l'ouverture du fichier de donnees :"<<filename<< endl;
        exit(-1);
        }
        // On lit le nombre d'echantillons
        double nbreEchantillon; 
        fscanf(pfile,"%le",&nbreEchantillon);
        nbSample=(long) nbreEchantillon;
        Transmi= NR.vector(1,nbSample);
        frequency=NR.vector(1,nbSample);
        
        double freq, filtretrans;
        for(int i=1; i<=nbSample; i++)
        {       fscanf(pfile,"%le%le",&freq,&filtretrans);
                frequency[i]=freq;
                Transmi[i]=filtretrans;
        }
        
        _numin=1.e99;
        _numax=0.;
        for(int i=1; i<=nbSample; i++)
        {       if(_numin>frequency[i]) _numin=frequency[i];
                if(_numax<frequency[i]) _numax=frequency[i];
        }
}

TabulatedFilter::~TabulatedFilter()     
{       NR.free_vector(Transmi,1,nbSample);
        NR.free_vector(frequency,1,nbSample);
}

double TabulatedFilter::transmission(double freq) const
{       
        if((freq<_numin)||(freq>_numax)) return 0.;
        
        NR.hunt(frequency,nbSample,freq,&LastfreqIndex);
        if((LastfreqIndex==0)||(LastfreqIndex==nbSample)) return 0.;
        else
        {       double Transinf=Transmi[LastfreqIndex];
                double Transsup=Transmi[LastfreqIndex+1];
                double slope=(Transsup-Transinf)/(frequency[LastfreqIndex+1]-frequency[LastfreqIndex]);
                double retVal=Transinf+slope*(freq-frequency[LastfreqIndex]);
                return retVal;
        }
}       
 
#include <fstream>

void TabulatedFilter::TraceFiltre(ostream& strOut, long nbPoint)
{       cout<<"On trace un TabulatedFilter"<<endl;
        strOut<<"Frequence(Hz)"<<'\t'<<"Transmission"<<endl;
        
        double pasFreq= (_numax-_numin)/nbPoint;
        double freq;
        for(long i=0; i<nbPoint; i++)
        {       freq=_numin+pasFreq*i;
                strOut<<freq<<'\t'<<transmission(freq)<<endl;
        }
        return;
}