|
5
|
1 /*
|
|
|
2 * SpanDSP - a series of DSP components for telephony
|
|
|
3 *
|
|
|
4 * awgn.c - An additive Gaussian white noise generator
|
|
|
5 *
|
|
|
6 * Written by Steve Underwood <steveu@coppice.org>
|
|
|
7 *
|
|
|
8 * Copyright (C) 2001 Steve Underwood
|
|
|
9 *
|
|
|
10 * All rights reserved.
|
|
|
11 *
|
|
|
12 * This program is free software; you can redistribute it and/or modify
|
|
|
13 * it under the terms of the GNU General Public License version 2, as
|
|
|
14 * published by the Free Software Foundation.
|
|
|
15 *
|
|
|
16 * This program is distributed in the hope that it will be useful,
|
|
|
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
19 * GNU General Public License for more details.
|
|
|
20 *
|
|
|
21 * You should have received a copy of the GNU General Public License
|
|
|
22 * along with this program; if not, write to the Free Software
|
|
|
23 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
|
24 *
|
|
|
25 * $Id: awgn.c,v 1.11 2006/11/19 14:07:24 steveu Exp $
|
|
|
26 */
|
|
|
27
|
|
|
28 /*! \file */
|
|
|
29
|
|
|
30 /* This code is based on some demonstration code in a research
|
|
|
31 paper somewhere. I can't track down where I got the original from,
|
|
|
32 so that due recognition can be given. The original had no explicit
|
|
|
33 copyright notice, and I hope nobody objects to its use here.
|
|
|
34
|
|
|
35 Having a reasonable Gaussian noise generator is pretty important for
|
|
|
36 telephony testing (in fact, pretty much any DSP testing), and this
|
|
|
37 one seems to have served me OK. Since the generation of Gaussian
|
|
|
38 noise is only for test purposes, and not a core system component,
|
|
|
39 I don't intend to worry excessively about copyright issues, unless
|
|
|
40 someone worries me.
|
|
|
41
|
|
|
42 The non-core nature of this code also explains why it is unlikely
|
|
|
43 to ever be optimised. */
|
|
|
44
|
|
|
45 #ifdef HAVE_CONFIG_H
|
|
|
46 #include <config.h>
|
|
|
47 #endif
|
|
|
48
|
|
|
49 #include <stdlib.h>
|
|
|
50 #include <inttypes.h>
|
|
|
51 #if defined(HAVE_TGMATH_H)
|
|
|
52 #include <tgmath.h>
|
|
|
53 #endif
|
|
|
54 #if defined(HAVE_MATH_H)
|
|
|
55 #include <math.h>
|
|
|
56 #endif
|
|
|
57
|
|
|
58 #include "spandsp/telephony.h"
|
|
|
59 #include "spandsp/dc_restore.h"
|
|
|
60 #include "spandsp/awgn.h"
|
|
|
61
|
|
|
62 /* Gaussian noise generator constants */
|
|
|
63 #define M1 259200
|
|
|
64 #define IA1 7141
|
|
|
65 #define IC1 54773
|
|
|
66 #define RM1 (1.0/M1)
|
|
|
67 #define M2 134456
|
|
|
68 #define IA2 8121
|
|
|
69 #define IC2 28411
|
|
|
70 #define RM2 (1.0/M2)
|
|
|
71 #define M3 243000
|
|
|
72 #define IA3 4561
|
|
|
73 #define IC3 51349
|
|
|
74
|
|
|
75 static double ran1(awgn_state_t *s)
|
|
|
76 {
|
|
|
77 double temp;
|
|
|
78 int j;
|
|
|
79
|
|
|
80 s->ix1 = (IA1*s->ix1 + IC1)%M1;
|
|
|
81 s->ix2 = (IA2*s->ix2 + IC2)%M2;
|
|
|
82 s->ix3 = (IA3*s->ix3 + IC3)%M3;
|
|
|
83 j = 1 + ((97*s->ix3)/M3);
|
|
|
84 if (j > 97 || j < 1)
|
|
|
85 {
|
|
|
86 /* Error */
|
|
|
87 return -1;
|
|
|
88 }
|
|
|
89 temp = s->r[j];
|
|
|
90 s->r[j] = (s->ix1 + s->ix2*RM2)*RM1;
|
|
|
91 return temp;
|
|
|
92 }
|
|
|
93 /*- End of function --------------------------------------------------------*/
|
|
|
94
|
|
|
95 void awgn_init_dbm0(awgn_state_t *s, int idum, float level)
|
|
|
96 {
|
|
|
97 awgn_init_dbov(s, idum, level - DBM0_MAX_POWER);
|
|
|
98 }
|
|
|
99 /*- End of function --------------------------------------------------------*/
|
|
|
100
|
|
|
101 void awgn_init_dbov(awgn_state_t *s, int idum, float level)
|
|
|
102 {
|
|
|
103 int j;
|
|
|
104
|
|
|
105 if (idum < 0)
|
|
|
106 idum = -idum;
|
|
|
107
|
|
|
108 s->rms = pow(10.0, level/20.0)*32768.0;
|
|
|
109
|
|
|
110 s->ix1 = (IC1 + idum)%M1;
|
|
|
111 s->ix1 = (IA1*s->ix1 + IC1)%M1;
|
|
|
112 s->ix2 = s->ix1%M2;
|
|
|
113 s->ix1 = (IA1*s->ix1 + IC1)%M1;
|
|
|
114 s->ix3 = s->ix1%M3;
|
|
|
115 s->r[0] = 0.0;
|
|
|
116 for (j = 1; j <= 97; j++)
|
|
|
117 {
|
|
|
118 s->ix1 = (IA1*s->ix1 + IC1)%M1;
|
|
|
119 s->ix2 = (IA2*s->ix2 + IC2)%M2;
|
|
|
120 s->r[j] = (s->ix1 + s->ix2*RM2)*RM1;
|
|
|
121 }
|
|
|
122 s->gset = 0.0;
|
|
|
123 s->iset = 0;
|
|
|
124 }
|
|
|
125 /*- End of function --------------------------------------------------------*/
|
|
|
126
|
|
|
127 int16_t awgn(awgn_state_t *s)
|
|
|
128 {
|
|
|
129 double fac;
|
|
|
130 double r;
|
|
|
131 double v1;
|
|
|
132 double v2;
|
|
|
133 double amp;
|
|
|
134
|
|
|
135 if (s->iset == 0)
|
|
|
136 {
|
|
|
137 do
|
|
|
138 {
|
|
|
139 v1 = 2.0*ran1(s) - 1.0;
|
|
|
140 v2 = 2.0*ran1(s) - 1.0;
|
|
|
141 r = v1*v1 + v2*v2;
|
|
|
142 }
|
|
|
143 while (r >= 1.0);
|
|
|
144 fac = sqrt(-2.0*log(r)/r);
|
|
|
145 s->gset = v1*fac;
|
|
|
146 s->iset = 1;
|
|
|
147 amp = v2*fac*s->rms;
|
|
|
148 }
|
|
|
149 else
|
|
|
150 {
|
|
|
151 s->iset = 0;
|
|
|
152 amp = s->gset*s->rms;
|
|
|
153 }
|
|
|
154 return fsaturate(amp);
|
|
|
155 }
|
|
|
156 /*- End of function --------------------------------------------------------*/
|
|
|
157 /*- End of file ------------------------------------------------------------*/
|
