Discussion :

[kchiheb]

comment ajouter un bruit gaussien a une fonction

[FR119492]

Salut!
J'ai retrouvé des sous-programmes Fortran permettant de générer un bruit gaussien. Il te suffira de les traduire.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
      SUBROUTINE GGNML  (DSEED,NR,R)
C
C   PURPOSE             - NORMAL OR GAUSSIAN RANDOM DEVIATE GENERATOR
C
C   USAGE               - CALL GGNML (DSEED,NR,R)
C
C   ARGUMENTS    DSEED  - INPUT/OUTPUT DOUBLE PRECISION VARIABLE
C                           ASSIGNED AN INTEGER VALUE IN THE
C                           EXCLUSIVE RANGE (1.D0, 2147483647.D0).
C                           DSEED IS REPLACED BY A NEW VALUE TO BE
C                           USED IN A SUBSEQUENT CALL.
C                NR     - INPUT NUMBER OF DEVIATES TO BE GENERATED.
C                R      - OUTPUT VECTOR OF LENGTH NR CONTAINING THE
C                           NORMAL (0,1) RANDOM NUMBERS.
C
      INTEGER            NR
      REAL               R(NR)
      DOUBLE PRECISION   DSEED
C
      INTEGER             IER
C
      CALL GGUBS(DSEED,NR,R)
C                                  TRANSFORMS EACH UNIFORM DEVIATE
      DO 5 I=1,NR
         CALL MDNRIS(R(I),R(I),IER)
    5 CONTINUE
      RETURN
      END

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
      SUBROUTINE GGUBS (DSEED,NR,R)
C
C   ARGUMENTS    DSEED  - INPUT/OUTPUT DOUBLE PRECISION VARIABLE
C                           ASSIGNED AN INTEGER VALUE IN THE
C                           EXCLUSIVE RANGE (1.D0, 2147483647.D0).
C                           DSEED IS REPLACED BY A NEW VALUE TO BE
C                           USED IN A SUBSEQUENT CALL.
C                NR     - INPUT NUMBER OF DEVIATES TO BE GENERATED.
C                R      - OUTPUT VECTOR OF LENGTH NR CONTAINING THE
C                           PSEUDO-RANDOM UNIFORM (0,1) DEVIATES
C
      INTEGER            NR
      REAL               R(NR)
      DOUBLE PRECISION   DSEED
      INTEGER            I
      DOUBLE PRECISION   D2P31M,D2P31
      DATA               D2P31M/2147483647.D0/
      DATA               D2P31/2147483648.D0/
C
      DO 5 I=1,NR
         DSEED = DMOD(16807.D0*DSEED,D2P31M)
    5 R(I) = DSEED / D2P31
      RETURN
      END

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
      SUBROUTINE MDNRIS (P,Y,IER)
C
C   PURPOSE             - INVERSE STANDARD NORMAL (GAUSSIAN)
C                           PROBABILITY DISTRIBUTION FUNCTION
C
C   ARGUMENTS    P      - INPUT VALUE IN THE EXCLUSIVE RANGE (0.0,1.0)
C                Y      - OUTPUT VALUE OF THE INVERSE NORMAL (0,1)
C                           PROBABILITY DISTRIBUTION FUNCTION
C                IER    - ERROR PARAMETER (OUTPUT)
C                         TERMINAL ERROR
C                           IER = 129 INDICATES P LIES OUTSIDE THE LEGAL
C                             RANGE. PLUS OR MINUS MACHINE INFINITY IS
C                             GIVEN AS THE RESULT (SIGN IS THE SIGN OF
C                             THE FUNCTION VALUE OF THE NEAREST LEGAL
C                             ARGUMENT).
C
C
      REAL               P,Y
      INTEGER            IER
      REAL               EPS,G0,G1,G2,G3,H0,H1,H2,A,W,WI,SN,SD
      REAL               SIGMA,SQRT2,X,XINF
      DATA               XINF/0.3402E+39/
      DATA               SQRT2/1.414214/
      DATA               EPS/0.9537E-06/
      DATA               G0/.1851159E-3/,G1/-.2028152E-2/
      DATA               G2/-.1498384/,G3/.1078639E-1/
      DATA               H0/.9952975E-1/,H1/.5211733/
      DATA               H2/-.6888301E-1/
C
      IER = 0
      IF (P .GT. 0.0 .AND. P .LT. 1.0) GO TO 5
      IER = 129
      SIGMA = SIGN(1.0,P)
      Y = SIGMA * XINF
      GO TO 9000
    5 IF(P.LE.EPS) GO TO 10
      X = 1.0 -(P + P)
      CALL MERFI (X,Y,IER)
      Y = -SQRT2 * Y
      GO TO 9005
C
   10 A = P+P
      W = SQRT(-ALOG(A+(A-A*A)))
C
      WI = 1./W
      SN = ((G3*WI+G2)*WI+G1)*WI
      SD = ((WI+H2)*WI+H1)*WI+H0
      Y = W + W*(G0+SN/SD)
      Y = -Y*SQRT2
      GO TO 9005
 9000 CONTINUE
      CALL UERTST(IER,6HMDNRIS)
 9005 RETURN
      END

Jean-Marc Blanc

[mr_samurai]

Salut,

Tu as regardé la fonction RANDN ? Elle génère un série suivant une loi normale centrée réduite:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
A = A + randn(size(A));

++