Mercurial > hg > audiostuff

comparison spandsp-0.0.3/spandsp-0.0.3/src/awgn.c @ 5:f762bf195c4b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import spandsp-0.0.3
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Fri, 25 Jun 2010 16:00:21 +0200
parents
children
comparison
equal deleted inserted replaced
4:26cd8f1ef0b1 5:f762bf195c4b
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 ------------------------------------------------------------*/

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