source: Sophya/trunk/Poubelle/DPC:FitsIOServer/Samba/bruit.cc@ 2995

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#ifdef __MWERKS__
   #include "mwerksmath.h"
   #include "unixmac.h"
#endif

#include "fmath.h"

#include "bruit.h"
#include "nbrandom.h"
// #include "rancern.h"
// #include "hbook.h"

//  Le code des classes   NoiseGenerator  RWalkNoise

//++ 
// Class        NoiseGenerator
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links        Childs
//
//    RWalkNoise OOFNoise EXPNoise MemNoise  SumNoise
//     
//-- 
//++
// Titre        Constructors
//--

/* --Methode-- */
//++
NoiseGenerator::NoiseGenerator(float sigma)
//
//--
{
if (sigma < 0.)  sigma = 1.;
mNCoups = 0;
mSigma = sigma;
//printf("-- NoiseGenerator::NoiseGenerator(%g) (Constructeur) ---\n", sigma);
}

//++
// Titre        Destructor
//--
/* --Methode-- */
//++
NoiseGenerator::~NoiseGenerator()
//
//--
{
    //printf("-- NoiseGenerator::~NoiseGenerator() (Destructeur) --- \n");
}


//++
//
// inline  unsigned long int  NoiseGenerator::NCoups() 
//--
//++
// Titre        Public Methods
//--

/* --Methode-- */
//++
float NoiseGenerator::Noise()
//
//--
{
mNCoups++;
return(NorRand()*mSigma);
}
//++
// Class        RWalkNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links        Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre        Constructor
//--

/* --Methode-- */
//++
RWalkNoise::RWalkNoise(float sigma)
    : NoiseGenerator(sigma)
//
//--
{
mState = 0.;
//printf("--  RWalkNoise::RWalkNoise(%g) (Constructeur) ---\n", sigma);
}


//++
// Titre        Destructor
//--
/* --Methode-- */
//++
RWalkNoise::~RWalkNoise()
//
//--
{
//printf("--  RWalkNoise::~RWalkNoise (Destructeur) ---\n");

}

//++
// Titre        Public Methods
//--
/* --Methode-- */
//++
float RWalkNoise::Noise()
//
//--
{
mState += NoiseGenerator::Noise();
return(mState);
}

//++
// Class        OOFNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links        Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre        Constructor
//--

/* --Methode-- */
//++
OOFNoise::OOFNoise(float sigma, int typacf, int mem, float tau)
    : NoiseGenerator(sigma)
//
//--
{
if (typacf != ACF_Exp)   typacf = ACF_Exp;
// Fonction d'autocorrelation donnant un bruit en 1/f
//                  
mTypACF = typacf;
if (mem < 1)  mem = 1;
mMemL = mem;
if (tau < 1.e-6)  tau = 1.e-6;
mTau = tau;
mState = new double[mem];
mDyn = new double[mem];
int i;
mDyn[0]=1.;
mState[0]=0.;
for(i=1; i<mem; i++) {
  mState[i] = 0.;   mDyn[i] = 1./sqrt((double)i);
}
//printf("--  OOFNoise::OOFNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
//       sigma, typacf, mem, tau);
}


//++
// Titre        Destructor
//--
/* --Methode-- */
//++
OOFNoise::~OOFNoise()
//
//--
{
delete[] mState;
delete[] mDyn;
// printf("--  OOFNoise::~OOFNoise (Destructeur) ---\n");
}


//++
// Titre        Public Methods
//--
/* --Methode-- */
//++
float OOFNoise::Noise()
//
//--
{
int i;
double rn;

for(i=(mMemL-1); i>0; i--)   mState[i] = mState[i-1];
mState[0] = NoiseGenerator::Noise();
rn = 0.;
for(i=0; i< mMemL; i++)  rn += mDyn[i]*mState[i];
return(rn);
}

/* --Methode-- */
//++
void OOFNoise::Print()
//
//--
{
int i,j;
printf("OOFNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
printf("OOFNoise::Print() Vecteur de Dynamique / Etat : \n");
for(i=0; i<mMemL-1; i+=2) 
  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
return;
}
//++
// Class        EXPNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links        Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre        Constructor
//--

/* --Methode-- */
//++
EXPNoise::EXPNoise(float sigma, int typacf, int mem, float tau)
    : NoiseGenerator(sigma)
//
//--
{
if (typacf != ACF_Exp)   typacf = ACF_Exp;
// Fonction d'autocorrelation exponentiel disponible uniquement
//                  Reza , Juillet 97
mTypACF = typacf;
if (mem < 1)  mem = 1;
mMemL = mem;
if (tau < 1.e-6)  tau = 1.e-6;
mTau = tau;
mState = new double[mem];
mDyn = new double[mem];
int i;
for(i=0; i<mem; i++) {
  mState[i] = 0.;   mDyn[i] = exp(-(double)i/tau);
}
//printf("--  EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
//       sigma, typacf, mem, tau);
}


//++
// Titre        Destructor
//--
/* --Methode-- */
//++
EXPNoise::~EXPNoise()
//
//--
{
delete[] mState;
delete[] mDyn;
//printf("--  EXPNoise::~EXPNoise (Destructeur) ---\n");
}


//++
// Titre        Public Methods
//--
/* --Methode-- */
//++
float EXPNoise::Noise()
//
//--
{
int i;
double rn;

for(i=(mMemL-1); i>0; i--)   mState[i] = mState[i-1];
mState[0] = NoiseGenerator::Noise();
rn = 0.;
for(i=0; i< mMemL; i++)  rn += mDyn[i]*mState[i];
return(rn);
}

/* --Methode-- */
//++
void EXPNoise::Print()
//
//--
{
int i,j;
printf("EXPNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
printf("EXPNoise::Print() Vecteur de Dynamique / Etat : \n");
for(i=0; i<mMemL-1; i+=2) 
  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   

return;
}

//++
// Class        MemNoise
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links        Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre        Constructor
//--



/* --Methode-- */
//++
MemNoise::MemNoise(float sigma, int mem, float tau, int ava)
    : NoiseGenerator(sqrtf(fabsf(tau))*sigma)
//
//--
{
/* on tire les instants des impulsions successives selon une
loi de distance exponentielle de duree tau
*/
const int place = 1000;
mTypACF = ACF_Exp;
if (mem < 1)  mem = 1;
mMemL = place;
mCoupHaut=mem;
if (tau < 1.)  tau = 1.;
mTau = tau;
mStPos = new float[place];
mStNeg = new float[place];
mTePos = new float[place];
mTeNeg = new float[place];
// printf("apres reservation\n");
int i;
for(i=0; i<place; i++)
{
mStNeg[i] = 0.; mStPos[i] = 0.;
}
// printf("apres initialisation tableaux \n");

mNappel=0;
mNtirage=1;
mNappLast=0;
mMemPos=1;
mMemNeg=1;
mduree=0.;
mTePos[0]=0.;
mTeNeg[0]=0.;
mTdernier=0.;

Avance(ava);
// printf("fin du createur MemNoise \n");
//printf("--  EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n", 
//       sigma, typacf, mem, tau);
}


//++
// Titre        Destructor
//--
/* --Methode-- */
//++
MemNoise::~MemNoise()
//
//--
{
delete[] mStPos;
delete[] mStNeg;
delete[] mTePos;
delete[] mTeNeg;
//printf("--  EXPNoise::~EXPNoise (Destructeur) ---\n");
}

//++
// Titre        Public Methods
//--
/* --Methode-- */
//++
float MemNoise::Noise()
//
//--
{
return(Avance(0));
}

/* --Methode-- */
//++
float MemNoise::Avance(long Asauter)
//
//--
{
int i,j;
float boum;
float rn;
float Somme,SbaryT,Poids,BaryT,Approx,Erreur;
int nReduit;
const float epsilon=.05;
// printf("entree dans methode Noise de MemNoise \n");
long mNsaut;
for (mNsaut = Asauter ; mNsaut >= 0 ; mNsaut--)
    {  // debut mNsaut
        // printf("mNsaut = %ld\n",mNsaut);
mNappel++;
mTdernier += 1.;
for(i=0; i<mMemPos; i++) mTePos[i] += 1.;
for(i=0; i<mMemNeg; i++) mTeNeg[i] += 1.;

if (mduree < mTdernier)
    {
        /* la duree depuis le dernier tirage est maintenant ecoulee
           on commence par remettre a jour le temps ecoule entre la derniere
           impulsion (celle qu'on va tirer maintenant) et l'instant present  */
           
           mTdernier = mTdernier - mduree;
        
        /* puis on tire une nouvelle impulsion de bruit */

           boum =  NoiseGenerator::Noise();

        /* on decale la memoire d'une case vers le passe 
           puis on met a jour la case [0] */

if(boum>0)
    {
        for(i=mMemPos; i>0; i--)
            {
                mStPos[i] = mStPos[i-1];
                mTePos[i] = mTePos[i-1];
            }
        mTePos[0] = mTdernier;
        mStPos[0] = boum;
        if (mMemPos<mMemL) mMemPos++;
/* a chaque tirage, on regroupe les impulsions precedentes */
for(i=mMemPos-2 ; i>0 ; i--)
    {
        if(mStPos[i]==0.)continue;
        //      printf("i,temps,bruit %d,%g,%g \n",i,mTePos[i],mStPos[i]);
        Somme = mStPos[i]/sqrtf(mTePos[i]);
            SbaryT = mStPos[i] * mTePos[i];
            Poids = mStPos[i];
            for(j=i-1 ; j>=0 ; j--)
                {
                    //          printf("  j,temps,bruit %d,%g,%g \n",j,mTePos[j],mStPos[j]);
                Somme += mStPos[j]/sqrtf(mTePos[j]);
                    Poids += mStPos[j];
                    SbaryT += mStPos[j] * mTePos[j];
                    BaryT = SbaryT / Poids;
                    Approx = Poids / sqrtf(BaryT);
                    Erreur = (Approx - Somme) / Somme;
                    //                              printf("i,j = %d %d Approx = %g \n",i,j,fabsf(Erreur));
                    //              printf("Temps,Bruit = %g %g \n",BaryT,Poids);
                    if(fabsf(Erreur)<epsilon)
                        {
                            mStPos[i]=Poids;
                            mTePos[i]=BaryT;
                            mStPos[j]=0.;
                        }
                    else break;
                }
    }
nReduit=0;
for(i=0; i<mMemPos-1 ; i++)
    {
        if(mStPos[i]!=0)
            {
                mStPos[nReduit]=mStPos[i];
                mTePos[nReduit]=mTePos[i];
                nReduit++;
            }
    }
mMemPos=nReduit+1;

    }
else
    {
        for(i=mMemNeg; i>0; i--)
            {
                mStNeg[i] = mStNeg[i-1];
                mTeNeg[i] = mTeNeg[i-1];
            }
        mTeNeg[0] = mTdernier;
        mStNeg[0] = boum;
        if (mMemNeg<mMemL) mMemNeg++;
    

/* idem pour negatifs */
for(i=mMemNeg-2 ; i>0 ; i--)
    {
        if(mStNeg[i]==0.)continue;
        //      printf("i,temps,bruit %d,%g,%g \n",i,mTeNeg[i],mStNeg[i]);
        Somme = mStNeg[i]/sqrtf(mTeNeg[i]);
            SbaryT = mStNeg[i] * mTeNeg[i];
            Poids = mStNeg[i];
            for(j=i-1 ; j>=0 ; j--)
                {
                    //          printf("  j,temps,bruit %d,%g,%g \n",j,mTeNeg[j],mStNeg[j]);
                Somme += mStNeg[j]/sqrtf(mTeNeg[j]);
                    Poids += mStNeg[j];
                    SbaryT += mStNeg[j] * mTeNeg[j];
                    BaryT = SbaryT / Poids;
                    Approx = Poids / sqrtf(BaryT);
                    Erreur = (Approx - Somme) / Somme;
                    //              printf("i,j = %d %d Approx = %g \n",i,j,Erreur);
                    //              printf("Temps,Bruit = %g %g \n",BaryT,Poids);
                    if(fabsf(Erreur)<epsilon)
                        {
                            mStNeg[i]=Poids;
                            mTeNeg[i]=BaryT;
                            mStNeg[j]=0.;
                        }
                    else break;
                }
    }
nReduit=0;
for(i=0; i<mMemNeg-1 ; i++)
    {
        if(mStNeg[i]!=0)
            {
                mStNeg[nReduit]=mStNeg[i];
                mTeNeg[nReduit]=mTeNeg[i];
                nReduit++;
            }
    }
mMemNeg=nReduit+1;
    }
/* on tire la duree a attendre avant la prochaine impulsion de bruit */
// $CHECK$ avec Francois - rndm() remplace par Reza 10/01/99
mduree = -mTau*log(frand01());  

// if(mNtirage<mMemL) mNtirage++;
mNtirage++;
/*
int idn=1;
float pds=1.;
float nul=0.;
hfill_(&idn,&mduree,&nul,&pds);
*/
// printf("Ntirage %d,mMemPos,Neg %d,%d, Suppress = %d \n",mNtirage,mMemPos,mMemNeg, Suppress);
    }
/* On calcule le bruit total */

//printf("mMemPos,Neg,%d %d\n",mMemPos,mMemNeg);
//for(i=0;i<mMemPos;i++) printf("mTePos %g \n",mTePos);
//for(i=0;i<mMemNeg;i++) printf("mTeNeg %g \n",mTeNeg);
if(mNsaut == 0)
    {
        // printf("coucou\n");
rn = 0.;
int indlim=0;
for(i=0; i < mMemPos-1; i++) 
{
if(mTePos[i] > .1)break; 
rn += mStPos[i];
indlim++;
}
// float ajout=0.;
for(i=indlim; i < mMemPos-1; i++)  rn += mStPos[i]/sqrtf(mTePos[i]);
indlim=0;
for(i=0; i < mMemNeg-1; i++) 
{
if(mTeNeg[i] > .1) break;
rn += mStNeg[i];
indlim++;
}
for(i=indlim; i < mMemNeg-1; i++) rn += mStNeg[i]/sqrtf(mTeNeg[i]);
}
}
return(rn);
}
/* --Methode-- */
//++
int MemNoise::Print()
//
//--
{
int i,j,rc=0;
printf("MemNoise::Print()  MemL=%d  Tau= %g \n", mMemL, mTau);
printf("mNtirage= %d , mMemPos= %d , mMemNeg= %d \n",mNtirage,mMemPos,mMemNeg);
// printf("MemNoise::Print() Vecteur de Dynamique / Etat : \n");
/* for(i=0; i<mMemL-1; i+=2) 
  printf("%d:   D=%g  S=%g   | %d:  D=%g  S=%g \n", i, 
         mDyn[i], mState[i], i+1, mDyn[i], mState[i]);   
         */
return(rc);
}


//++
// Class        SumNoise
//    Generateur de bruit blanc + 1/f avec fknee
//
// include      bruit.h math.h fmath.h nbrandom.h
//
//--
//++
//
// Links        Parents
//
//    NoiseGenerator
//     
//-- 
//++
// Titre        Constructor
//--
// Generateur de bruit blanc + 1/f avec fknee


/* --Methode-- */
//++
SumNoise::SumNoise(float fknee, float sig)
    : NoiseGenerator(sig)
//
//--
{
// Reza 27/01/98 :
// Je calcule les sigma du MemNoise d'apres la formule de Francois 
mn = new MemNoise(sig*fknee);  
}


//++
// Titre        Destructor
//--
/* --Methode-- */
//++
SumNoise::~SumNoise()
//
//--
{
delete mn;
}

//++
// Titre        Public Methods
//--
/* --Methode-- */
//++
float SumNoise::Noise()
//
//--
{
return(NoiseGenerator::Noise() + mn->Noise());
}