1 /*
2 * SpanDSP - a series of DSP components for telephony
3 *
4 * makecss.c - Create the composite source signal (CSS) for G.168 testing.
5 *
6 * Written by Steve Underwood <steveu@coppice.org>
7 *
8 * Copyright (C) 2003 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
26 /*! \page makecss_page CSS construction for G.168 testing
27 \section makecss_page_sec_1 What does it do?
28 ???.
29
30 \section makecss_page_sec_2 How does it work?
31 ???.
32 */
33
34 #if defined(HAVE_CONFIG_H)
35 #include "config.h"
36 #endif
37
38 #include <stdlib.h>
39 #include <unistd.h>
40 #include <string.h>
41 #include <time.h>
42 #include <stdio.h>
43 #include <fcntl.h>
44 #include <sndfile.h>
45 #if defined(HAVE_FFTW3_H)
46 #include <fftw3.h>
47 #else
48 #include <fftw.h>
49 #endif
50
51 #include "spandsp.h"
52 #include "spandsp/g168models.h"
53
54 #if !defined(NULL)
55 #define NULL (void *) 0
56 #endif
57
58 #define FAST_SAMPLE_RATE 44100.0
59
60 #define C1_VOICED_SAMPLES 2144 /* 48.62ms at 44100 samples/second => 2144.142 */
61 #define C1_NOISE_SAMPLES 8820 /* 200ms at 44100 samples/second => 8820.0 */
62 #define C1_SILENCE_SAMPLES 4471 /* 101.38ms at 44100 samples/second => 4470.858 */
63
64 #define C3_VOICED_SAMPLES 3206 /* 72.69ms at 44100 samples/second => 3205.629 */
65 #define C3_NOISE_SAMPLES 8820 /* 200ms at 44100 samples/second => 8820.0 */
66 #define C3_SILENCE_SAMPLES 5614 /* 127.31ms at 44100 samples/second => 5614.371 */
67
scaling(double f,double start,double end,double start_gain,double end_gain)68 static double scaling(double f, double start, double end, double start_gain, double end_gain)
69 {
70 double scale;
71
72 scale = start_gain + (f - start)*(end_gain - start_gain)/(end - start);
73 return scale;
74 }
75 /*- End of function --------------------------------------------------------*/
76
peak(const int16_t amp[],int len)77 static double peak(const int16_t amp[], int len)
78 {
79 int16_t peak;
80 int i;
81
82 peak = 0;
83 for (i = 0; i < len; i++)
84 {
85 if (abs(amp[i]) > peak)
86 peak = abs(amp[i]);
87 }
88 return peak/32767.0;
89 }
90 /*- End of function --------------------------------------------------------*/
91
rms(const int16_t amp[],int len)92 static double rms(const int16_t amp[], int len)
93 {
94 double ms;
95 int i;
96
97 ms = 0.0;
98 for (i = 0; i < len; i++)
99 ms += amp[i]*amp[i];
100 return sqrt(ms/len)/32767.0;
101 }
102 /*- End of function --------------------------------------------------------*/
103
rms_to_dbm0(double rms)104 static double rms_to_dbm0(double rms)
105 {
106 return 20.0*log10(rms) + DBM0_MAX_POWER;
107 }
108 /*- End of function --------------------------------------------------------*/
109
rms_to_db(double rms)110 static double rms_to_db(double rms)
111 {
112 return 20.0*log10(rms);
113 }
114 /*- End of function --------------------------------------------------------*/
115
main(int argc,char * argv[])116 int main(int argc, char *argv[])
117 {
118 #if defined(HAVE_FFTW3_H)
119 double in[8192][2];
120 double out[8192][2];
121 #else
122 fftw_complex in[8192];
123 fftw_complex out[8192];
124 #endif
125 fftw_plan p;
126 int16_t voiced_sound[8192];
127 int16_t noise_sound[8830];
128 int16_t silence_sound[8192];
129 int i;
130 int voiced_length;
131 int randy;
132 double f;
133 double pk;
134 double ms;
135 double scale;
136 SNDFILE *filehandle;
137 SF_INFO info;
138 awgn_state_t *noise_source;
139
140 memset(&info, 0, sizeof(info));
141 info.frames = 0;
142 info.samplerate = FAST_SAMPLE_RATE;
143 info.channels = 1;
144 info.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
145 info.sections = 1;
146 info.seekable = 1;
147 if ((filehandle = sf_open("sound_c1.wav", SFM_WRITE, &info)) == NULL)
148 {
149 fprintf(stderr, " Failed to open result file\n");
150 exit(2);
151 }
152
153 printf("Generate C1\n");
154 /* The sequence is
155 48.62ms of voiced sound from table C.1 of G.168
156 200.0ms of pseudo-noise
157 101.38ms of silence
158 The above is then repeated phase inverted.
159
160 The voice comes straight from C.1, repeated enough times to
161 fill out the 48.62ms period - i.e. 16 copies of the sequence.
162
163 The pseudo noise section is random numbers filtered by the spectral
164 pattern in Figure C.2 */
165
166 /* The set of C1 voice samples is ready for use in the output file. */
167 voiced_length = sizeof(css_c1)/sizeof(css_c1[0]);
168 for (i = 0; i < voiced_length; i++)
169 voiced_sound[i] = css_c1[i];
170 pk = peak(voiced_sound, voiced_length);
171 ms = rms(voiced_sound, voiced_length);
172 printf("Voiced level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));
173
174 #if defined(HAVE_FFTW3_H)
175 p = fftw_plan_dft_1d(8192, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
176 #else
177 p = fftw_create_plan(8192, FFTW_BACKWARD, FFTW_ESTIMATE);
178 #endif
179 for (i = 0; i < 8192; i++)
180 {
181 #if defined(HAVE_FFTW3_H)
182 in[i][0] = 0.0;
183 in[i][1] = 0.0;
184 #else
185 in[i].re = 0.0;
186 in[i].im = 0.0;
187 #endif
188 }
189 for (i = 1; i <= 3715; i++)
190 {
191 f = FAST_SAMPLE_RATE*i/8192.0;
192
193 #if 1
194 if (f < 50.0)
195 scale = -60.0;
196 else if (f < 100.0)
197 scale = scaling(f, 50.0, 100.0, -25.8, -12.8);
198 else if (f < 200.0)
199 scale = scaling(f, 100.0, 200.0, -12.8, 17.4);
200 else if (f < 215.0)
201 scale = scaling(f, 200.0, 215.0, 17.4, 17.8);
202 else if (f < 500.0)
203 scale = scaling(f, 215.0, 500.0, 17.8, 12.2);
204 else if (f < 1000.0)
205 scale = scaling(f, 500.0, 1000.0, 12.2, 7.2);
206 else if (f < 2850.0)
207 scale = scaling(f, 1000.0, 2850.0, 7.2, 0.0);
208 else if (f < 3600.0)
209 scale = scaling(f, 2850.0, 3600.0, 0.0, -2.0);
210 else if (f < 3660.0)
211 scale = scaling(f, 3600.0, 3660.0, -2.0, -20.0);
212 else if (f < 3680.0)
213 scale = scaling(f, 3600.0, 3680.0, -20.0, -30.0);
214 else
215 scale = -60.0;
216 #else
217 scale = 0.0;
218 #endif
219 randy = ((rand() >> 10) & 0x1) ? 1.0 : -1.0;
220 #if defined(HAVE_FFTW3_H)
221 in[i][0] = randy*pow(10.0, scale/20.0)*35.0;
222 in[8192 - i][0] = -in[i][0];
223 #else
224 in[i].re = randy*pow(10.0, scale/20.0)*35.0;
225 in[8192 - i].re = -in[i].re;
226 #endif
227 }
228 #if defined(HAVE_FFTW3_H)
229 fftw_execute(p);
230 #else
231 fftw_one(p, in, out);
232 #endif
233 for (i = 0; i < 8192; i++)
234 {
235 #if defined(HAVE_FFTW3_H)
236 noise_sound[i] = out[i][1];
237 #else
238 noise_sound[i] = out[i].im;
239 #endif
240 }
241 pk = peak(noise_sound, 8192);
242 ms = rms(noise_sound, 8192);
243 printf("Filtered noise level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));
244
245 for (i = 0; i < 8192; i++)
246 silence_sound[i] = 0.0;
247
248 for (i = 0; i < 16; i++)
249 sf_writef_short(filehandle, voiced_sound, voiced_length);
250 printf("%d samples of voice\n", 16*voiced_length);
251 sf_writef_short(filehandle, noise_sound, 8192);
252 sf_writef_short(filehandle, noise_sound, C1_NOISE_SAMPLES - 8192);
253 printf("%d samples of noise\n", C1_NOISE_SAMPLES);
254 sf_writef_short(filehandle, silence_sound, C1_SILENCE_SAMPLES);
255 printf("%d samples of silence\n", C1_SILENCE_SAMPLES);
256
257 /* Now phase invert the C1 set of samples. */
258 voiced_length = sizeof(css_c1)/sizeof(css_c1[0]);
259 for (i = 0; i < voiced_length; i++)
260 voiced_sound[i] = -css_c1[i];
261 for (i = 0; i < 8192; i++)
262 noise_sound[i] = -noise_sound[i];
263
264 for (i = 0; i < 16; i++)
265 sf_writef_short(filehandle, voiced_sound, voiced_length);
266 printf("%d samples of voice\n", 16*voiced_length);
267 sf_writef_short(filehandle, noise_sound, 8192);
268 sf_writef_short(filehandle, noise_sound, C1_NOISE_SAMPLES - 8192);
269 printf("%d samples of noise\n", C1_NOISE_SAMPLES);
270 sf_writef_short(filehandle, silence_sound, C1_SILENCE_SAMPLES);
271 printf("%d samples of silence\n", C1_SILENCE_SAMPLES);
272
273 if (sf_close(filehandle))
274 {
275 fprintf(stderr, " Cannot close speech file '%s'\n", "sound_c1.wav");
276 exit(2);
277 }
278
279 memset(&info, 0, sizeof(info));
280 info.frames = 0;
281 info.samplerate = FAST_SAMPLE_RATE;
282 info.channels = 1;
283 info.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;
284 info.sections = 1;
285 info.seekable = 1;
286 if ((filehandle = sf_open("sound_c3.wav", SFM_WRITE, &info)) == NULL)
287 {
288 fprintf(stderr, " Failed to open result file\n");
289 exit(2);
290 }
291
292 printf("Generate C3\n");
293 /* The sequence is
294 72.69ms of voiced sound from table C.3 of G.168
295 200.0ms of pseudo-noise
296 127.31ms of silence
297 The above is then repeated phase inverted.
298
299 The voice comes straight from C.3, repeated enough times to
300 fill out the 72.69ms period - i.e. 14 copies of the sequence.
301
302 The pseudo noise section is AWGN filtered by the spectral
303 pattern in Figure C.2. Since AWGN has the quality of being its
304 own Fourier transform, we can use an approach like the one above
305 for the C1 signal, using AWGN samples instead of randomly alternating
306 ones and zeros. */
307
308 /* Take the supplied set of C3 voice samples. */
309 voiced_length = (sizeof(css_c3)/sizeof(css_c3[0]));
310 for (i = 0; i < voiced_length; i++)
311 voiced_sound[i] = css_c3[i];
312 pk = peak(voiced_sound, voiced_length);
313 ms = rms(voiced_sound, voiced_length);
314 printf("Voiced level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));
315
316 noise_source = awgn_init_dbm0(NULL, 7162534, rms_to_dbm0(ms));
317 for (i = 0; i < 8192; i++)
318 noise_sound[i] = awgn(noise_source);
319 pk = peak(noise_sound, 8192);
320 ms = rms(noise_sound, 8192);
321 printf("Unfiltered noise level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));
322
323 /* Now filter them */
324 #if defined(HAVE_FFTW3_H)
325 p = fftw_plan_dft_1d(8192, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
326 #else
327 p = fftw_create_plan(8192, FFTW_BACKWARD, FFTW_ESTIMATE);
328 #endif
329 for (i = 0; i < 8192; i++)
330 {
331 #if defined(HAVE_FFTW3_H)
332 in[i][0] = 0.0;
333 in[i][1] = 0.0;
334 #else
335 in[i].re = 0.0;
336 in[i].im = 0.0;
337 #endif
338 }
339 for (i = 1; i <= 3715; i++)
340 {
341 f = FAST_SAMPLE_RATE*i/8192.0;
342
343 #if 1
344 if (f < 50.0)
345 scale = -60.0;
346 else if (f < 100.0)
347 scale = scaling(f, 50.0, 100.0, -25.8, -12.8);
348 else if (f < 200.0)
349 scale = scaling(f, 100.0, 200.0, -12.8, 17.4);
350 else if (f < 215.0)
351 scale = scaling(f, 200.0, 215.0, 17.4, 17.8);
352 else if (f < 500.0)
353 scale = scaling(f, 215.0, 500.0, 17.8, 12.2);
354 else if (f < 1000.0)
355 scale = scaling(f, 500.0, 1000.0, 12.2, 7.2);
356 else if (f < 2850.0)
357 scale = scaling(f, 1000.0, 2850.0, 7.2, 0.0);
358 else if (f < 3600.0)
359 scale = scaling(f, 2850.0, 3600.0, 0.0, -2.0);
360 else if (f < 3660.0)
361 scale = scaling(f, 3600.0, 3660.0, -2.0, -20.0);
362 else if (f < 3680.0)
363 scale = scaling(f, 3600.0, 3680.0, -20.0, -30.0);
364 else
365 scale = -60.0;
366 #else
367 scale = 0.0;
368 #endif
369 #if defined(HAVE_FFTW3_H)
370 in[i][0] = noise_sound[i]*pow(10.0, scale/20.0)*0.0106;
371 in[8192 - i][0] = -in[i][0];
372 #else
373 in[i].re = noise_sound[i]*pow(10.0, scale/20.0)*0.0106;
374 in[8192 - i].re = -in[i].re;
375 #endif
376 }
377 #if defined(HAVE_FFTW3_H)
378 fftw_execute(p);
379 #else
380 fftw_one(p, in, out);
381 #endif
382 for (i = 0; i < 8192; i++)
383 {
384 #if defined(HAVE_FFTW3_H)
385 noise_sound[i] = out[i][1];
386 #else
387 noise_sound[i] = out[i].im;
388 #endif
389 }
390 pk = peak(noise_sound, 8192);
391 ms = rms(noise_sound, 8192);
392 printf("Filtered noise level = %.2fdB, crest factor = %.2fdB\n", rms_to_dbm0(ms), rms_to_db(pk/ms));
393
394 for (i = 0; i < 14; i++)
395 sf_writef_short(filehandle, voiced_sound, voiced_length);
396 printf("%d samples of voice\n", 14*voiced_length);
397 sf_writef_short(filehandle, noise_sound, 8192);
398 sf_writef_short(filehandle, noise_sound, C3_NOISE_SAMPLES - 8192);
399 printf("%d samples of noise\n", C3_NOISE_SAMPLES);
400 sf_writef_short(filehandle, silence_sound, C3_SILENCE_SAMPLES);
401 printf("%d samples of silence\n", C3_SILENCE_SAMPLES);
402
403 /* Now phase invert the C3 set of samples. */
404 voiced_length = (sizeof(css_c3)/sizeof(css_c3[0]));
405 for (i = 0; i < voiced_length; i++)
406 voiced_sound[i] = -css_c3[i];
407 for (i = 0; i < 8192; i++)
408 noise_sound[i] = -noise_sound[i];
409
410 for (i = 0; i < 14; i++)
411 sf_writef_short(filehandle, voiced_sound, voiced_length);
412 printf("%d samples of voice\n", 14*voiced_length);
413 sf_writef_short(filehandle, noise_sound, 8192);
414 sf_writef_short(filehandle, noise_sound, C3_NOISE_SAMPLES - 8192);
415 printf("%d samples of noise\n", C3_NOISE_SAMPLES);
416 sf_writef_short(filehandle, silence_sound, C3_SILENCE_SAMPLES);
417 printf("%d samples of silence\n", C3_SILENCE_SAMPLES);
418
419 if (sf_close(filehandle))
420 {
421 fprintf(stderr, " Cannot close speech file '%s'\n", "sound_c3.wav");
422 exit(2);
423 }
424
425 fftw_destroy_plan(p);
426 return 0;
427 }
428 /*- End of function --------------------------------------------------------*/
429 /*- End of file ------------------------------------------------------------*/
430