comparison spandsp-0.0.3/spandsp-0.0.3/tests/noise_tests.c @ 5:f762bf195c4b

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 * noise_tests.c
5 *
6 * Written by Steve Underwood <steveu@coppice.org>
7 *
8 * Copyright (C) 2005 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: noise_tests.c,v 1.9 2006/11/19 14:07:27 steveu Exp $
26 */
27
28 /*! \page noise_tests_page Noise generator tests
29 \section noise_tests_page_sec_1 What does it do?
30 */
31
32 #ifdef HAVE_CONFIG_H
33 #include "config.h"
34 #endif
35
36 #include <stdio.h>
37 #include <inttypes.h>
38 #include <stdlib.h>
39 #include <string.h>
40 #if defined(HAVE_TGMATH_H)
41 #include <tgmath.h>
42 #endif
43 #if defined(HAVE_MATH_H)
44 #include <math.h>
45 #endif
46 #include <audiofile.h>
47 #include <tiffio.h>
48
49 #include "spandsp.h"
50
51 #if !defined(M_PI)
52 # define M_PI 3.14159265358979323846 /* pi */
53 #endif
54
55 #define OUT_FILE_NAME "noise.wav"
56
57 /* Some simple sanity tests for the noise generation routines */
58
59 int main (int argc, char *argv[])
60 {
61 int i;
62 int j;
63 int level;
64 int clip_high;
65 int clip_low;
66 int total_samples;
67 int seed = 1234567;
68 int outframes;
69 int16_t value;
70 double total;
71 double x;
72 double p;
73 double o;
74 int bins[65536];
75 int16_t amp[1024];
76 noise_state_t noise_source;
77 AFfilehandle outhandle;
78 AFfilesetup filesetup;
79
80 if ((filesetup = afNewFileSetup()) == AF_NULL_FILESETUP)
81 {
82 fprintf(stderr, " Failed to create file setup\n");
83 exit(2);
84 }
85 afInitSampleFormat(filesetup, AF_DEFAULT_TRACK, AF_SAMPFMT_TWOSCOMP, 16);
86 afInitRate(filesetup, AF_DEFAULT_TRACK, (float) SAMPLE_RATE);
87 afInitFileFormat(filesetup, AF_FILE_WAVE);
88 afInitChannels(filesetup, AF_DEFAULT_TRACK, 1);
89 if ((outhandle = afOpenFile(OUT_FILE_NAME, "w", filesetup)) == AF_NULL_FILEHANDLE)
90 {
91 fprintf(stderr, " Cannot create wave file '%s'\n", OUT_FILE_NAME);
92 exit(2);
93 }
94
95 /* Generate AWGN at several RMS levels between -50dBOv and 0dBOv. Noise is
96 generated for a large number of samples (1,000,000), and the RMS value
97 of the noise is calculated along the way. If the resulting level is
98 close to the requested RMS level, at least the scaling of the noise
99 should be Ok. At high levels some clipping may distort the result a
100 little. */
101 printf("Testing with quality 7 AWGN\n");
102 for (level = -50; level <= 0; level += 5)
103 {
104 clip_high = 0;
105 clip_low = 0;
106 total = 0.0;
107 noise_init_dbov(&noise_source, seed, (float) level, NOISE_CLASS_AWGN, 7);
108 total_samples = 1000000;
109 for (i = 0; i < total_samples; i++)
110 {
111 value = noise(&noise_source);
112 if (value == 32767)
113 clip_high++;
114 else if (value == -32768)
115 clip_low++;
116 total += ((double) value)*((double) value);
117 }
118 printf ("RMS = %.3f (expected %d) %.2f%% error [clipped samples %d+%d]\n",
119 log10(sqrt(total/total_samples)/32768.0)*20.0,
120 level,
121 100.0*(1.0 - sqrt(total/total_samples)/(pow(10.0, level/20.0)*32768.0)),
122 clip_low,
123 clip_high);
124 if (level < -5 && fabs(log10(sqrt(total/total_samples)/32768.0)*20.0 - level) > 0.2)
125 {
126 printf("Test failed\n");
127 exit(2);
128 }
129 }
130
131 printf("Testing with quality 20 AWGN\n");
132 for (level = -50; level <= 0; level += 5)
133 {
134 clip_high = 0;
135 clip_low = 0;
136 total = 0.0;
137 noise_init_dbov(&noise_source, seed, (float) level, NOISE_CLASS_AWGN, 20);
138 total_samples = 1000000;
139 for (i = 0; i < total_samples; i++)
140 {
141 value = noise(&noise_source);
142 if (value == 32767)
143 clip_high++;
144 else if (value == -32768)
145 clip_low++;
146 total += ((double) value)*((double) value);
147 }
148 printf ("RMS = %.3f (expected %d) %.2f%% error [clipped samples %d+%d]\n",
149 log10(sqrt(total/total_samples)/32768.0)*20.0,
150 level,
151 100.0*(1.0 - sqrt(total/total_samples)/(pow(10.0, level/20.0)*32768.0)),
152 clip_low,
153 clip_high);
154 if (level < -5 && fabs(log10(sqrt(total/total_samples)/32768.0)*20.0 - level) > 0.2)
155 {
156 printf("Test failed\n");
157 exit(2);
158 }
159 }
160
161 /* Now look at the statistical spread of the results, by collecting data in
162 bins from a large number of samples. Use a fairly high noise level, but
163 low enough to avoid significant clipping. Use the Gaussian model to
164 predict the real probability, and present the results for graphing. */
165 memset(bins, 0, sizeof(bins));
166 clip_high = 0;
167 clip_low = 0;
168 level = -15;
169 noise_init_dbov(&noise_source, seed, (float) level, NOISE_CLASS_AWGN, 7);
170 total_samples = 10000000;
171 for (i = 0; i < total_samples; i++)
172 {
173 value = noise(&noise_source);
174 if (value == 32767)
175 clip_high++;
176 else if (value == -32768)
177 clip_low++;
178 bins[value + 32768]++;
179 }
180 /* Find the RMS power level to expect */
181 o = pow(10.0, level/20.0)*(32768.0*0.70711);
182 for (i = 0; i < 65536 - 10; i++)
183 {
184 x = i - 32768;
185 /* Find the real probability for this bin */
186 p = (1.0/(o*sqrt(2.0*M_PI)))*exp(-(x*x)/(2.0*o*o));
187 /* Now do a little smoothing on the real data to get a reasonably
188 steady answer */
189 x = 0;
190 for (j = 0; j < 10; j++)
191 x += bins[i + j];
192 x /= 10.0;
193 x /= total_samples;
194 /* Now send it out for graphing. */
195 if (p > 0.0000001)
196 printf("%6d %.7f %.7f\n", i - 32768, x, p);
197 }
198
199 printf("Generating AWGN at -15dBOv to file\n");
200 for (j = 0; j < 50; j++)
201 {
202 for (i = 0; i < 1024; i++)
203 amp[i] = noise(&noise_source);
204 outframes = afWriteFrames(outhandle,
205 AF_DEFAULT_TRACK,
206 amp,
207 1024);
208 if (outframes != 1024)
209 {
210 fprintf(stderr, " Error writing wave file\n");
211 exit(2);
212 }
213 }
214
215 /* Generate AWGN at several RMS levels between -50dBm and 0dBm. Noise is
216 generated for a large number of samples (1,000,000), and the RMS value
217 of the noise is calculated along the way. If the resulting level is
218 close to the requested RMS level, at least the scaling of the noise
219 should be Ok. At high levels some clipping may distort the result a
220 little. */
221 printf("Testing with quality 7 Hoth noise.\n");
222 for (level = -50; level <= 0; level += 5)
223 {
224 clip_high = 0;
225 clip_low = 0;
226 total = 0.0;
227 noise_init_dbov(&noise_source, seed, (float) level, NOISE_CLASS_HOTH, 7);
228 total_samples = 1000000;
229 for (i = 0; i < total_samples; i++)
230 {
231 value = noise(&noise_source);
232 if (value == 32767)
233 clip_high++;
234 else if (value == -32768)
235 clip_low++;
236 total += ((double) value)*((double) value);
237 }
238 printf ("RMS = %.3f (expected %d) %.2f%% error [clipped samples %d+%d]\n",
239 log10(sqrt(total/total_samples)/32768.0)*20.0,
240 level,
241 100.0*(1.0 - sqrt(total/total_samples)/(pow(10.0, level/20.0)*32768.0)),
242 clip_low,
243 clip_high);
244 if (level < -5 && fabs(log10(sqrt(total/total_samples)/32768.0)*20.0 - level) > 0.2)
245 {
246 printf("Test failed\n");
247 exit(2);
248 }
249 }
250
251 printf("Generating Hoth noise at -15dBOv to file\n");
252 level = -15;
253 noise_init_dbov(&noise_source, seed, (float) level, NOISE_CLASS_HOTH, 7);
254 for (j = 0; j < 50; j++)
255 {
256 for (i = 0; i < 1024; i++)
257 amp[i] = noise(&noise_source);
258 outframes = afWriteFrames(outhandle,
259 AF_DEFAULT_TRACK,
260 amp,
261 1024);
262 if (outframes != 1024)
263 {
264 fprintf(stderr, " Error writing wave file\n");
265 exit(2);
266 }
267 }
268
269 if (afCloseFile(outhandle))
270 {
271 fprintf(stderr, " Cannot close wave file '%s'\n", OUT_FILE_NAME);
272 exit(2);
273 }
274 afFreeFileSetup(filesetup);
275
276 printf("Tests passed.\n");
277 return 0;
278 }
279 /*- End of function --------------------------------------------------------*/
280 /*- End of file ------------------------------------------------------------*/

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