source: Sophya/trunk/SophyaLib/Samba/bruit.cc@ 4035

Last change on this file since 4035 was 3572, checked in by cmv, 17 years ago

char* -> const char* pour regler les problemes de deprecated string const... + comparaison unsigned signed + suppression EVOL_PLANCK rz+cmv 07/02/2009

File size: 12.5 KB
Line 
1#include <stdlib.h>
2#include <stdio.h>
3#include <math.h>
4
5#include "sopnamsp.h"
6#include "fmath.h"
7
8#include "bruit.h"
9#include "srandgen.h"
10// #include "rancern.h"
11// #include "hbook.h"
12#ifdef Darwin
13#include <limits.h>
14#endif
15
16// Le code des classes NoiseGenerator RWalkNoise
17
18//++
19// Class NoiseGenerator
20//
21// include bruit.h math.h fmath.h srandgen.h
22//
23//--
24//++
25//
26// Links Childs
27//
28// RWalkNoise OOFNoise EXPNoise MemNoise SumNoise
29//
30//--
31//++
32// Titre Constructors
33//--
34
35// Definition de fonctions pour float - a changer en utilisant cmath
36// Reza Juillet 2002
37inline float _myfabsf(float x) { return ((float)(fabs((double)x)) ); }
38inline float _mysqrtf(float x) { return ((float)(sqrt((double)x)) ); }
39
40/* --Methode-- */
41//++
42NoiseGenerator::NoiseGenerator(float sigma)
43//
44//--
45{
46if (sigma < 0.) sigma = 1.;
47mNCoups = 0;
48mSigma = sigma;
49//printf("-- NoiseGenerator::NoiseGenerator(%g) (Constructeur) ---\n", sigma);
50}
51
52//++
53// Titre Destructor
54//--
55/* --Methode-- */
56//++
57NoiseGenerator::~NoiseGenerator()
58//
59//--
60{
61 //printf("-- NoiseGenerator::~NoiseGenerator() (Destructeur) --- \n");
62}
63
64
65//++
66//
67// inline unsigned long int NoiseGenerator::NCoups()
68//--
69//++
70// Titre Public Methods
71//--
72
73/* --Methode-- */
74//++
75float NoiseGenerator::Noise()
76//
77//--
78{
79mNCoups++;
80return(NorRand()*mSigma);
81}
82//++
83// Class RWalkNoise
84//
85// include bruit.h math.h fmath.h srandgen.h
86//
87//--
88//++
89//
90// Links Parents
91//
92// NoiseGenerator
93//
94//--
95//++
96// Titre Constructor
97//--
98
99/* --Methode-- */
100//++
101RWalkNoise::RWalkNoise(float sigma)
102 : NoiseGenerator(sigma)
103//
104//--
105{
106mState = 0.;
107//printf("-- RWalkNoise::RWalkNoise(%g) (Constructeur) ---\n", sigma);
108}
109
110
111//++
112// Titre Destructor
113//--
114/* --Methode-- */
115//++
116RWalkNoise::~RWalkNoise()
117//
118//--
119{
120//printf("-- RWalkNoise::~RWalkNoise (Destructeur) ---\n");
121
122}
123
124//++
125// Titre Public Methods
126//--
127/* --Methode-- */
128//++
129float RWalkNoise::Noise()
130//
131//--
132{
133mState += NoiseGenerator::Noise();
134return(mState);
135}
136
137//++
138// Class OOFNoise
139//
140// include bruit.h math.h fmath.h srandgen.h
141//
142//--
143//++
144//
145// Links Parents
146//
147// NoiseGenerator
148//
149//--
150//++
151// Titre Constructor
152//--
153
154/* --Methode-- */
155//++
156OOFNoise::OOFNoise(float sigma, int typacf, int mem, float tau)
157 : NoiseGenerator(sigma)
158//
159//--
160{
161if (typacf != ACF_Exp) typacf = ACF_Exp;
162// Fonction d'autocorrelation donnant un bruit en 1/f
163//
164mTypACF = typacf;
165if (mem < 1) mem = 1;
166mMemL = mem;
167if (tau < 1.e-6) tau = 1.e-6;
168mTau = tau;
169mState = new double[mem];
170mDyn = new double[mem];
171int i;
172mDyn[0]=1.;
173mState[0]=0.;
174for(i=1; i<mem; i++) {
175 mState[i] = 0.; mDyn[i] = 1./sqrt((double)i);
176}
177//printf("-- OOFNoise::OOFNoise(%g,%d,%d,%g) (Constructeur) ---\n",
178// sigma, typacf, mem, tau);
179}
180
181
182//++
183// Titre Destructor
184//--
185/* --Methode-- */
186//++
187OOFNoise::~OOFNoise()
188//
189//--
190{
191delete[] mState;
192delete[] mDyn;
193// printf("-- OOFNoise::~OOFNoise (Destructeur) ---\n");
194}
195
196
197//++
198// Titre Public Methods
199//--
200/* --Methode-- */
201//++
202float OOFNoise::Noise()
203//
204//--
205{
206int i;
207double rn;
208
209for(i=(mMemL-1); i>0; i--) mState[i] = mState[i-1];
210mState[0] = NoiseGenerator::Noise();
211rn = 0.;
212for(i=0; i< mMemL; i++) rn += mDyn[i]*mState[i];
213return(rn);
214}
215
216/* --Methode-- */
217//++
218void OOFNoise::Print()
219//
220//--
221{
222printf("OOFNoise::Print() MemL=%d Tau= %g \n", mMemL, mTau);
223printf("OOFNoise::Print() Vecteur de Dynamique / Etat : \n");
224for(int i=0; i<mMemL-1; i+=2)
225 printf("%d: D=%g S=%g | %d: D=%g S=%g \n", i,
226 mDyn[i], mState[i], i+1, mDyn[i], mState[i]);
227return;
228}
229//++
230// Class EXPNoise
231//
232// include bruit.h math.h fmath.h srandgen.h
233//
234//--
235//++
236//
237// Links Parents
238//
239// NoiseGenerator
240//
241//--
242//++
243// Titre Constructor
244//--
245
246/* --Methode-- */
247//++
248EXPNoise::EXPNoise(float sigma, int typacf, int mem, float tau)
249 : NoiseGenerator(sigma)
250//
251//--
252{
253if (typacf != ACF_Exp) typacf = ACF_Exp;
254// Fonction d'autocorrelation exponentiel disponible uniquement
255// Reza , Juillet 97
256mTypACF = typacf;
257if (mem < 1) mem = 1;
258mMemL = mem;
259if (tau < 1.e-6) tau = 1.e-6;
260mTau = tau;
261mState = new double[mem];
262mDyn = new double[mem];
263int i;
264for(i=0; i<mem; i++) {
265 mState[i] = 0.; mDyn[i] = exp(-(double)i/tau);
266}
267//printf("-- EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n",
268// sigma, typacf, mem, tau);
269}
270
271
272//++
273// Titre Destructor
274//--
275/* --Methode-- */
276//++
277EXPNoise::~EXPNoise()
278//
279//--
280{
281delete[] mState;
282delete[] mDyn;
283//printf("-- EXPNoise::~EXPNoise (Destructeur) ---\n");
284}
285
286
287//++
288// Titre Public Methods
289//--
290/* --Methode-- */
291//++
292float EXPNoise::Noise()
293//
294//--
295{
296int i;
297double rn;
298
299for(i=(mMemL-1); i>0; i--) mState[i] = mState[i-1];
300mState[0] = NoiseGenerator::Noise();
301rn = 0.;
302for(i=0; i< mMemL; i++) rn += mDyn[i]*mState[i];
303return(rn);
304}
305
306/* --Methode-- */
307//++
308void EXPNoise::Print()
309//
310//--
311{
312printf("EXPNoise::Print() MemL=%ld Tau= %g \n", mMemL, mTau);
313printf("EXPNoise::Print() Vecteur de Dynamique / Etat : \n");
314for(long int i=0; i<mMemL-1; i+=2)
315 printf("%ld: D=%g S=%g | %ld: D=%g S=%g \n", i,
316 mDyn[i], mState[i], i+1, mDyn[i], mState[i]);
317
318return;
319}
320
321//++
322// Class MemNoise
323//
324// include bruit.h math.h fmath.h srandgen.h
325//
326//--
327//++
328//
329// Links Parents
330//
331// NoiseGenerator
332//
333//--
334//++
335// Titre Constructor
336//--
337
338
339// ---------------- $CHECK$ Reza 1/12/99 ------------
340// ----- MAJ MemNoise et SumNoise / version F. Couchot (~couchot/CoSa/Samba/bruit.cc )
341// -----------------------------------------------------------
342
343/* --Methode-- */
344//++
345MemNoise::MemNoise(float sigma, int mem, float tau, int ava)
346 : NoiseGenerator(_mysqrtf(_myfabsf(tau))*sigma)
347//
348//--
349{
350/* on tire les instants des impulsions successives selon une
351loi de distance exponentielle de duree tau
352*/
353const int place = 1000;
354mTypACF = ACF_Exp;
355if (mem < 1) mem = 1;
356mMemL = place;
357mCoupHaut=mem;
358if (tau < 1.) tau = 1.;
359mTau = tau;
360mStPos = new float[place];
361mStNeg = new float[place];
362mTePos = new float[place];
363mTeNeg = new float[place];
364// printf("apres reservation\n");
365int i;
366for(i=0; i<place; i++)
367{
368mStNeg[i] = 0.; mStPos[i] = 0.;
369}
370// printf("apres initialisation tableaux \n");
371
372mNappel=0;
373mNtirage=1;
374mNappLast=0;
375mMemPos=1;
376mMemNeg=1;
377mduree=0.;
378mTePos[0]=0.;
379mTeNeg[0]=0.;
380mTdernier=0.;
381
382Avance(ava);
383// printf("fin du createur MemNoise \n");
384//printf("-- EXPNoise::EXPNoise(%g,%d,%d,%g) (Constructeur) ---\n",
385// sigma, typacf, mem, tau);
386}
387
388
389//++
390// Titre Destructor
391//--
392/* --Methode-- */
393//++
394MemNoise::~MemNoise()
395//
396//--
397{
398delete[] mStPos;
399delete[] mStNeg;
400delete[] mTePos;
401delete[] mTeNeg;
402//printf("-- EXPNoise::~EXPNoise (Destructeur) ---\n");
403}
404
405//++
406// Titre Public Methods
407//--
408/* --Methode-- */
409//++
410float MemNoise::Noise()
411//
412//--
413{
414return(Avance(0));
415}
416
417/* --Methode-- */
418//++
419float MemNoise::Avance(long Asauter)
420//
421//--
422{
423int i,j;
424float boum;
425float rn = 0;
426float Somme,SbaryT,Poids,BaryT,Approx,Erreur;
427int nReduit;
428const float epsilon=.05;
429// printf("entree dans methode Noise de MemNoise \n");
430long mNsaut;
431for (mNsaut = Asauter ; mNsaut >= 0 ; mNsaut--)
432 { // debut mNsaut
433 // printf("mNsaut = %ld\n",mNsaut);
434mNappel++;
435mTdernier += 1.;
436for(i=0; i<mMemPos; i++) mTePos[i] += 1.;
437for(i=0; i<mMemNeg; i++) mTeNeg[i] += 1.;
438
439if (mduree < mTdernier)
440 {
441 /* la duree depuis le dernier tirage est maintenant ecoulee
442 on commence par remettre a jour le temps ecoule entre la derniere
443 impulsion (celle qu'on va tirer maintenant) et l'instant present */
444
445 mTdernier = mTdernier - mduree;
446
447 /* puis on tire une nouvelle impulsion de bruit */
448
449 boum = NoiseGenerator::Noise();
450
451 /* on decale la memoire d'une case vers le passe
452 puis on met a jour la case [0] */
453
454if(boum>0)
455 {
456 for(i=mMemPos; i>0; i--)
457 {
458 mStPos[i] = mStPos[i-1];
459 mTePos[i] = mTePos[i-1];
460 }
461 mTePos[0] = mTdernier;
462 mStPos[0] = boum;
463 if (mMemPos<mMemL) mMemPos++;
464/* a chaque tirage, on regroupe les impulsions precedentes */
465for(i=mMemPos-2 ; i>0 ; i--)
466 {
467 if(mStPos[i]==0.)continue;
468 // printf("i,temps,bruit %d,%g,%g \n",i,mTePos[i],mStPos[i]);
469 Somme = mStPos[i]/_mysqrtf(mTePos[i]);
470 SbaryT = mStPos[i] * mTePos[i];
471 Poids = mStPos[i];
472 for(j=i-1 ; j>=0 ; j--)
473 {
474 // printf(" j,temps,bruit %d,%g,%g \n",j,mTePos[j],mStPos[j]);
475 Somme += mStPos[j]/_mysqrtf(mTePos[j]);
476 Poids += mStPos[j];
477 SbaryT += mStPos[j] * mTePos[j];
478 BaryT = SbaryT / Poids;
479 Approx = Poids / _mysqrtf(BaryT);
480 Erreur = (Approx - Somme) / Somme;
481 // printf("i,j = %d %d Approx = %g \n",i,j,_myfabsf(Erreur));
482 // printf("Temps,Bruit = %g %g \n",BaryT,Poids);
483 if(_myfabsf(Erreur)<epsilon)
484 {
485 mStPos[i]=Poids;
486 mTePos[i]=BaryT;
487 mStPos[j]=0.;
488 }
489 else break;
490 }
491 }
492nReduit=0;
493for(i=0; i<mMemPos-1 ; i++)
494 {
495 if(mStPos[i]!=0)
496 {
497 mStPos[nReduit]=mStPos[i];
498 mTePos[nReduit]=mTePos[i];
499 nReduit++;
500 }
501 }
502mMemPos=nReduit+1;
503
504 }
505else
506 {
507 for(i=mMemNeg; i>0; i--)
508 {
509 mStNeg[i] = mStNeg[i-1];
510 mTeNeg[i] = mTeNeg[i-1];
511 }
512 mTeNeg[0] = mTdernier;
513 mStNeg[0] = boum;
514 if (mMemNeg<mMemL) mMemNeg++;
515
516
517/* idem pour negatifs */
518for(i=mMemNeg-2 ; i>0 ; i--)
519 {
520 if(mStNeg[i]==0.)continue;
521 // printf("i,temps,bruit %d,%g,%g \n",i,mTeNeg[i],mStNeg[i]);
522 Somme = mStNeg[i]/_mysqrtf(mTeNeg[i]);
523 SbaryT = mStNeg[i] * mTeNeg[i];
524 Poids = mStNeg[i];
525 for(j=i-1 ; j>=0 ; j--)
526 {
527 // printf(" j,temps,bruit %d,%g,%g \n",j,mTeNeg[j],mStNeg[j]);
528 Somme += mStNeg[j]/_mysqrtf(mTeNeg[j]);
529 Poids += mStNeg[j];
530 SbaryT += mStNeg[j] * mTeNeg[j];
531 BaryT = SbaryT / Poids;
532 Approx = Poids / _mysqrtf(BaryT);
533 Erreur = (Approx - Somme) / Somme;
534 // printf("i,j = %d %d Approx = %g \n",i,j,Erreur);
535 // printf("Temps,Bruit = %g %g \n",BaryT,Poids);
536 if(_myfabsf(Erreur)<epsilon)
537 {
538 mStNeg[i]=Poids;
539 mTeNeg[i]=BaryT;
540 mStNeg[j]=0.;
541 }
542 else break;
543 }
544 }
545nReduit=0;
546for(i=0; i<mMemNeg-1 ; i++)
547 {
548 if(mStNeg[i]!=0)
549 {
550 mStNeg[nReduit]=mStNeg[i];
551 mTeNeg[nReduit]=mTeNeg[i];
552 nReduit++;
553 }
554 }
555mMemNeg=nReduit+1;
556 }
557/* on tire la duree a attendre avant la prochaine impulsion de bruit */
558// $CHECK$ avec Francois - rndm() remplace par Reza 10/01/99
559// passage en drand le 1/12/99
560mduree = -mTau*log(drand01());
561
562// if(mNtirage<mMemL) mNtirage++;
563mNtirage++;
564/*
565int idn=1;
566float pds=1.;
567float nul=0.;
568hfill_(&idn,&mduree,&nul,&pds);
569*/
570// printf("Ntirage %d,mMemPos,Neg %d,%d, Suppress = %d \n",mNtirage,mMemPos,mMemNeg, Suppress);
571 }
572/* On calcule le bruit total */
573
574//printf("mMemPos,Neg,%d %d\n",mMemPos,mMemNeg);
575//for(i=0;i<mMemPos;i++) printf("mTePos %g \n",mTePos);
576//for(i=0;i<mMemNeg;i++) printf("mTeNeg %g \n",mTeNeg);
577if(mNsaut == 0)
578 {
579 // printf("coucou\n");
580rn = 0.;
581int indlim=0;
582for(i=0; i < mMemPos-1; i++)
583{
584if(mTePos[i] > 1.) break; // Modif F.C et D.Y. 03/2000
585rn += mStPos[i];
586indlim++;
587}
588// float ajout=0.;
589for(i=indlim; i < mMemPos-1; i++) rn += mStPos[i]/_mysqrtf(mTePos[i]);
590indlim=0;
591for(i=0; i < mMemNeg-1; i++)
592{
593if(mTeNeg[i] > 1.) break; // Modif F.C et D.Y. 03/2000
594rn += mStNeg[i];
595indlim++;
596}
597for(i=indlim; i < mMemNeg-1; i++) rn += mStNeg[i]/_mysqrtf(mTeNeg[i]);
598}
599}
600return(rn);
601}
602/* --Methode-- */
603//++
604int MemNoise::Print()
605//
606//--
607{
608int rc=0;
609printf("MemNoise::Print() MemL=%d Tau= %g \n", mMemL, mTau);
610printf("mNtirage= %d , mMemPos= %d , mMemNeg= %d \n",mNtirage,mMemPos,mMemNeg);
611// printf("MemNoise::Print() Vecteur de Dynamique / Etat : \n");
612/* for(int i=0; i<mMemL-1; i+=2)
613 printf("%d: D=%g S=%g | %d: D=%g S=%g \n", i,
614 mDyn[i], mState[i], i+1, mDyn[i], mState[i]);
615 */
616return(rc);
617}
618
619//++
620// Class SumNoise
621//
622// include bruit.h math.h fmath.h srandgen.h
623//
624// Generateur de bruit blanc + 1/f avec fknee
625
626//--
627//++
628//
629// Links Parents
630//
631// NoiseGenerator
632//
633//--
634//++
635// Titre Constructor
636//--
637// Generateur de bruit blanc + 1/f avec fknee
638
639
640/* --Methode-- */
641//++
642SumNoise::SumNoise(float fknee, float sig, int mem, float tau)
643 : NoiseGenerator(sig)
644//
645//--
646{
647// Reza 27/01/98 :
648// Je calcule les sigma du MemNoise d'apres la formule de Francois
649// mn = new MemNoise(sig*fknee);
650// mn = new MemNoise(1.4142*sig*_mysqrtf(fknee),1,10.,1000000); // modif appel FC 13 mai 98.
651mn = new MemNoise(1.4142*sig*sqrt((double)fknee),mem,tau,12000); // FC 13 mai 98. -> Reza 1/12/99
652}
653
654
655//++
656// Titre Destructor
657//--
658/* --Methode-- */
659//++
660SumNoise::~SumNoise()
661//
662//--
663{
664delete mn;
665}
666
667//++
668// Titre Public Methods
669//--
670/* --Methode-- */
671//++
672float SumNoise::Noise()
673//
674//--
675{
676return(NoiseGenerator::Noise() + mn->Noise());
677}
Note: See TracBrowser for help on using the repository browser.