Mercurial > hg > audiostuff

comparison spandsp-0.0.6pre17/src/awgn.c @ 4:26cd8f1ef0b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import spandsp-0.0.6pre17
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Fri, 25 Jun 2010 15:50:58 +0200
parents
children
comparison
equal deleted inserted replaced
3:c6c5a16ce2f2 4:26cd8f1ef0b1
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 Lesser General Public License version 2.1,
14 * as 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 Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License 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.22 2009/02/10 13:06:46 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 #if defined(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 #include "floating_fudge.h"
58
59 #include "spandsp/telephony.h"
60 #include "spandsp/fast_convert.h"
61 #include "spandsp/saturated.h"
62 #include "spandsp/awgn.h"
63
64 #include "spandsp/private/awgn.h"
65
66 /* Gaussian noise generator constants */
67 #define M1 259200
68 #define IA1 7141
69 #define IC1 54773
70 #define RM1 (1.0/M1)
71 #define M2 134456
72 #define IA2 8121
73 #define IC2 28411
74 #define RM2 (1.0/M2)
75 #define M3 243000
76 #define IA3 4561
77 #define IC3 51349
78
79 static double ran1(awgn_state_t *s)
80 {
81 double temp;
82 int j;
83
84 s->ix1 = (IA1*s->ix1 + IC1)%M1;
85 s->ix2 = (IA2*s->ix2 + IC2)%M2;
86 s->ix3 = (IA3*s->ix3 + IC3)%M3;
87 j = 1 + ((97*s->ix3)/M3);
88 if (j > 97 || j < 1)
89 {
90 /* Error */
91 return -1;
92 }
93 temp = s->r[j];
94 s->r[j] = (s->ix1 + s->ix2*RM2)*RM1;
95 return temp;
96 }
97 /*- End of function --------------------------------------------------------*/
98
99 SPAN_DECLARE(awgn_state_t *) awgn_init_dbov(awgn_state_t *s, int idum, float level)
100 {
101 int j;
102
103 if (s == NULL)
104 {
105 if ((s = (awgn_state_t *) malloc(sizeof(*s))) == NULL)
106 return NULL;
107 }
108 if (idum < 0)
109 idum = -idum;
110
111 s->rms = pow(10.0, level/20.0)*32768.0;
112
113 s->ix1 = (IC1 + idum)%M1;
114 s->ix1 = (IA1*s->ix1 + IC1)%M1;
115 s->ix2 = s->ix1%M2;
116 s->ix1 = (IA1*s->ix1 + IC1)%M1;
117 s->ix3 = s->ix1%M3;
118 s->r[0] = 0.0;
119 for (j = 1; j <= 97; j++)
120 {
121 s->ix1 = (IA1*s->ix1 + IC1)%M1;
122 s->ix2 = (IA2*s->ix2 + IC2)%M2;
123 s->r[j] = (s->ix1 + s->ix2*RM2)*RM1;
124 }
125 s->gset = 0.0;
126 s->iset = 0;
127 return s;
128 }
129 /*- End of function --------------------------------------------------------*/
130
131 SPAN_DECLARE(awgn_state_t *) awgn_init_dbm0(awgn_state_t *s, int idum, float level)
132 {
133 return awgn_init_dbov(s, idum, level - DBM0_MAX_POWER);
134 }
135 /*- End of function --------------------------------------------------------*/
136
137 SPAN_DECLARE(int) awgn_release(awgn_state_t *s)
138 {
139 return 0;
140 }
141 /*- End of function --------------------------------------------------------*/
142
143 SPAN_DECLARE(int) awgn_free(awgn_state_t *s)
144 {
145 free(s);
146 return 0;
147 }
148 /*- End of function --------------------------------------------------------*/
149
150 SPAN_DECLARE(int16_t) awgn(awgn_state_t *s)
151 {
152 double fac;
153 double r;
154 double v1;
155 double v2;
156 double amp;
157
158 if (s->iset == 0)
159 {
160 do
161 {
162 v1 = 2.0*ran1(s) - 1.0;
163 v2 = 2.0*ran1(s) - 1.0;
164 r = v1*v1 + v2*v2;
165 }
166 while (r >= 1.0);
167 fac = sqrt(-2.0*log(r)/r);
168 s->gset = v1*fac;
169 s->iset = 1;
170 amp = v2*fac*s->rms;
171 }
172 else
173 {
174 s->iset = 0;
175 amp = s->gset*s->rms;
176 }
177 return fsaturate(amp);
178 }
179 /*- End of function --------------------------------------------------------*/
180 /*- End of file ------------------------------------------------------------*/

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.