From e523020eb811b67bed175994a7f94c832652e157 Mon Sep 17 00:00:00 2001 From: Rick Sayre Date: Sat, 30 Dec 2023 19:11:48 -0800 Subject: [PATCH] speaker-test: Add bandwidth-limited pink noise at -18.5dB AES FS Based Closes: https://github.com/alsa-project/alsa-utils/pull/251 Signed-off-by: Rick Sayre Signed-off-by: Jaroslav Kysela --- speaker-test/Makefile.am | 4 +- speaker-test/speaker-test.1 | 9 +- speaker-test/speaker-test.c | 73 ++++++++++++---- speaker-test/st2095.c | 164 ++++++++++++++++++++++++++++++++++++ speaker-test/st2095.h | 41 +++++++++ 5 files changed, 272 insertions(+), 19 deletions(-) create mode 100644 speaker-test/st2095.c create mode 100644 speaker-test/st2095.h diff --git a/speaker-test/Makefile.am b/speaker-test/Makefile.am index d53dbd6..b21919b 100644 --- a/speaker-test/Makefile.am +++ b/speaker-test/Makefile.am @@ -3,7 +3,7 @@ SUBDIRS= samples LDADD = $(LIBINTL) -lm bin_PROGRAMS = speaker-test -speaker_test_SOURCES = speaker-test.c pink.c +speaker_test_SOURCES = speaker-test.c pink.c st2095.c man_MANS = speaker-test.1 -EXTRA_DIST = readme.txt speaker-test.1 pink.h +EXTRA_DIST = readme.txt speaker-test.1 pink.h st2095.h diff --git a/speaker-test/speaker-test.1 b/speaker-test/speaker-test.1 index add6b21..f375e19 100644 --- a/speaker-test/speaker-test.1 +++ b/speaker-test/speaker-test.1 @@ -94,11 +94,17 @@ Use number of periods. The default value is 4. stream of \fIRATE\fP Hz .TP -\fB\-t\fP | \fB\-\-test\fP \fBpink\fP|\fBsine\fP|\fBwav\fP +\fB\-t\fP | \fB\-\-test\fP \fBpink\fP|\fBst2095\fP|\fBsine\fP|\fBwav\fP \fB\-t pink\fP means use pink noise (default). Pink noise is perceptually uniform noise -- that is, it sounds like every frequency at once. If you can hear any tone it may indicate resonances in your speaker system or room. +\fB\-t st2095\fP means use bandlimited pink noise at -18.5dB AES FS, generated according to SMPTE ST-2095:1-2015. +In addition to speaker localization it may be used for system calibration, for example 85dB for thater drivers, with an extra +10dB for subwoofers. +Per the spec, it is intended "to be used in calibrating the sound pressure level and +electroacoustic response of a cinema B-chain system." +Note that sampling rates less than 48KHz are outside the scope of the spec, and an attempt will be made to construct a reduced rate filter. + \fB\-t sine\fP means to use sine wave. \fB\-t wav\fP means to play WAV files, either pre-defined files or given via \fB\-w\fP option. @@ -180,3 +186,4 @@ Play in the order of front\-right and front-left from the front PCM The speaker\-test program was written by James Courtier-Dutton. Pink noise support was added by Nathan Hurst. Further extensions by Takashi Iwai. +SMPTE ST-2095:1 band-limited pink noise added by Rick Sayre. diff --git a/speaker-test/speaker-test.c b/speaker-test/speaker-test.c index a13e737..1121401 100644 --- a/speaker-test/speaker-test.c +++ b/speaker-test/speaker-test.c @@ -25,8 +25,13 @@ * Some cleanup from Daniel Caujolle-Bert * Pink noise option added Nathan Hurst, * based on generator by Phil Burk (pink.c) + * ST-2095 noise option added Rick Sayre, + * based on generator specified by SMPTE ST-2095:1-2015 + * Also switched to stable harmonic oscillator for sine * * Changelog: + * 0.0.9 Added support for ST-2095 band-limited pink noise output, switched to harmonic oscillator for sine + * Changelog: * 0.0.8 Added support for pink noise output. * Changelog: * 0.0.7 Added support for more than 6 channels. @@ -55,6 +60,7 @@ #include #include #include "pink.h" +#include "st2095.h" #include "gettext.h" #include "version.h" #include "os_compat.h" @@ -71,6 +77,7 @@ enum { TEST_PINK_NOISE = 1, TEST_SINE, TEST_WAV, + TEST_ST2095_NOISE, TEST_PATTERN, }; @@ -373,16 +380,17 @@ static void do_generate(uint8_t *frames, int channel, int count, * Sine generator */ typedef struct { - double phase; - double max_phase; - double step; + double a; + double s; + double c; } sine_t; static void init_sine(sine_t *sine) { - sine->phase = 0; - sine->max_phase = 1.0 / freq; - sine->step = 1.0 / (double)rate; + // symplectic integration for fast, stable harmonic oscillator + sine->a = 2.0*M_PI * freq / rate; + sine->c = 1.0; + sine->s = 0.0; } static value_t generate_sine(void *arg) @@ -390,13 +398,13 @@ static value_t generate_sine(void *arg) sine_t *sine = arg; value_t res; - res.f = sin((sine->phase * 2 * M_PI) / sine->max_phase - M_PI); - res.f *= generator_scale; + res.f = sine->s * generator_scale; if (format != SND_PCM_FORMAT_FLOAT_LE) res.i = res.f * INT32_MAX; - sine->phase += sine->step; - if (sine->phase >= sine->max_phase) - sine->phase -= sine->max_phase; + + // update the oscillator + sine->c -= sine->a * sine->s; + sine->s += sine->a * sine->c; return res; } @@ -414,6 +422,20 @@ static value_t generate_pink_noise(void *arg) return res; } +/* Band-Limited Pink Noise, per SMPTE ST 2095-1 + * beyond speaker localization, this can be used for setting loudness to standard + */ +static value_t generate_st2095_noise(void *arg) +{ + st2095_noise_t *st2095 = arg; + value_t res; + + res.f = generate_st2095_noise_sample(st2095); + if (format != SND_PCM_FORMAT_FLOAT_LE) + res.i = res.f * INT32_MAX; + return res; +} + /* * useful for tests */ @@ -853,10 +875,14 @@ static int write_buffer(snd_pcm_t *handle, uint8_t *ptr, int cptr) static int pattern; static sine_t sine; static pink_noise_t pink; +static st2095_noise_t st2095; static void init_loop(void) { switch (test_type) { + case TEST_ST2095_NOISE: + initialize_st2095_noise(&st2095, rate); + break; case TEST_PINK_NOISE: initialize_pink_noise(&pink, 16); break; @@ -901,7 +927,12 @@ static int write_loop(snd_pcm_t *handle, int channel, int periods, uint8_t *fram do_generate(frames, channel, period_size, generate_pink_noise, &pink); else if (test_type == TEST_PATTERN) do_generate(frames, channel, period_size, generate_pattern, &pattern); - else + else if (test_type == TEST_ST2095_NOISE) { + reset_st2095_noise_measurement(&st2095); + do_generate(frames, channel, period_size, generate_st2095_noise, &st2095); + printf(_("\tSMPTE ST-2095 noise batch was %2.2fdB RMS\n"), + compute_st2095_noise_measurement(&st2095, period_size)); + } else do_generate(frames, channel, period_size, generate_sine, &sine); if ((err = write_buffer(handle, frames, period_size)) < 0) @@ -953,7 +984,7 @@ static void help(void) "-b,--buffer ring buffer size in us\n" "-p,--period period size in us\n" "-P,--nperiods number of periods\n" - "-t,--test pink=use pink noise, sine=use sine wave, wav=WAV file\n" + "-t,--test pink=use pink noise, sine=use sine wave, st2095=use SMPTE ST-2095 noise, wav=WAV file\n" "-l,--nloops specify number of loops to test, 0 = infinite\n" "-s,--speaker single speaker test. Values 1=Left, 2=right, etc\n" "-w,--wavfile Use the given WAV file as a test sound\n" @@ -1082,9 +1113,16 @@ int main(int argc, char *argv[]) { case 't': if (*optarg == 'p') test_type = TEST_PINK_NOISE; - else if (*optarg == 's') - test_type = TEST_SINE; - else if (*optarg == 'w') + else if (*optarg == 's') { + if (optarg[1] == 'i') + test_type = TEST_SINE; + else if (optarg[1] == 't') + test_type = TEST_ST2095_NOISE; + else { + fprintf(stderr, _("Invalid test type %s\n"), optarg); + exit(1); + } + } else if (*optarg == 'w') test_type = TEST_WAV; else if (*optarg == 't') test_type = TEST_PATTERN; @@ -1160,6 +1198,9 @@ int main(int argc, char *argv[]) { printf(_("Playback device is %s\n"), device); printf(_("Stream parameters are %iHz, %s, %i channels\n"), rate, snd_pcm_format_name(format), channels); switch (test_type) { + case TEST_ST2095_NOISE: + printf(_("Using SMPTE ST-2095 -18.5dB AES FS band-limited pink noise\n")); + break; case TEST_PINK_NOISE: printf(_("Using 16 octaves of pink noise\n")); break; diff --git a/speaker-test/st2095.c b/speaker-test/st2095.c new file mode 100644 index 0000000..f036074 --- /dev/null +++ b/speaker-test/st2095.c @@ -0,0 +1,164 @@ +/* + st2095.c + + Generate Bandlimited Pink Noise (-18.5dB AES FS) + Using the SMPTE ST 2095:1-2015 standard + + Based on pseudo-code from the above SMPTE standard, which bore the credit + "Revised 2015-01-04 by Calvert Dayton" + + Copyleft 2023 Rick Sayre - No rights reserved. +*/ + +#include "aconfig.h" +#include +#include +#include "st2095.h" + +/************************************************************/ + + +void reset_st2095_noise_measurement( st2095_noise_t *st2095 ) { + st2095->accum = 0.; +} + +float compute_st2095_noise_measurement( st2095_noise_t *st2095, int period ) { + return(10. * log10f(st2095->accum / (float)period) + 3.01); +} + +void initialize_st2095_noise( st2095_noise_t *st2095, int sample_rate) { + // Periodicity in samples must be a power of two, <= 2^31 + // Typical values are 524288, 1048576, 2097152 or 4194304 + if (sample_rate > 48000) { + // Special case LCG step for 1024K samples @ 88.2K or 96k + st2095->samplesPerPeriod = 1048576; + st2095->randStep = 163841; + } else { + st2095->samplesPerPeriod = 524288; + st2095->randStep = 52737; + } + + // set up LCG PRNG + st2095->randMax = st2095->samplesPerPeriod - 1; + st2095->seed = 0; + st2095->scaleFactor = 2.0 / (float)st2095->randMax; + + st2095->maxAmp = powf(10.0, ST2095_MAX_PEAK / 20.0); + + // Calculate omegaT for matched Z transform highpass filters + st2095->w0t = 2.0 * M_PI * ST2095_HPFC / (float)sample_rate; + + // Limit LpFc <= Nyquist (actually lower, based on 48 vs 22.4 KHz spec cutoff) + // The spec says the filter begins at 22.4KHz, if we ask for a Nyquist-impossible + // sampling rate, compute something with the same relationship + st2095->LpFc = ST2095_LPFC; + float rateratio = 48000. / ST2095_LPFC; + if (st2095->LpFc > sample_rate/rateratio) + st2095->LpFc = sample_rate/rateratio; + + // Calculate k and k^2 for bilinear transform lowpass filters + st2095->k = tanf(( 2.0 * M_PI * st2095->LpFc / (float)sample_rate ) / 2.0); + st2095->k2 = st2095->k * st2095->k; + + // Calculate biquad coefficients for bandpass filter components + st2095->hp1_a1 = -2.0 * expf(-0.3826835 * st2095->w0t) * cosf(0.9238795 * st2095->w0t); + st2095->hp1_a2 = expf(2.0 * -0.3826835 * st2095->w0t); + st2095->hp1_b0 = (1.0 - st2095->hp1_a1 + st2095->hp1_a2) / 4.0; + st2095->hp1_b1 = -2.0 * st2095->hp1_b0; + st2095->hp1_b2 = st2095->hp1_b0; + + st2095->hp2_a1 = -2.0 * expf(-0.9238795 * st2095->w0t) * cosf(0.3826835 * st2095->w0t); + st2095->hp2_a2 = expf(2.0 * -0.9238795 * st2095->w0t); + st2095->hp2_b0 = (1.0 - st2095->hp2_a1 + st2095->hp2_a2) / 4.0; + st2095->hp2_b1 = -2.0 * st2095->hp2_b0; + st2095->hp2_b2 = st2095->hp2_b0; + + st2095->lp1_a1 = (2.0 * (st2095->k2 - 1.0)) / (st2095->k2 + (st2095->k / 1.306563) + 1.0); + st2095->lp1_a2 = (st2095->k2 - (st2095->k / 1.306563) + 1.0) / (st2095->k2 + (st2095->k / 1.306563) + 1.0); + st2095->lp1_b0 = st2095->k2 / (st2095->k2 + (st2095->k / 1.306563) + 1.0); + st2095->lp1_b1 = 2.0 * st2095->lp1_b0; + st2095->lp1_b2 = st2095->lp1_b0; + + st2095->lp2_a1 = (2.0 * (st2095->k2 - 1.0)) / (st2095->k2 + (st2095->k / 0.541196) + 1.0); + st2095->lp2_a2 = (st2095->k2 - (st2095->k / 0.541196) + 1.0) / (st2095->k2 + (st2095->k / 0.541196) + 1.0); + st2095->lp2_b0 = st2095->k2 / (st2095->k2 + (st2095->k / 0.541196) + 1.0); + st2095->lp2_b1 = 2.0 * st2095->lp2_b0; + st2095->lp2_b2 = st2095->lp2_b0; + + // initialize delay lines for bandpass filter + st2095->hp1w1 = 0.0; + st2095->hp1w2 = 0.0; + st2095->hp2w1 = 0.0; + st2095->hp2w2 = 0.0; + st2095->lp1w1 = 0.0; + st2095->lp1w2 = 0.0; + st2095->lp2w1 = 0.0; + st2095->lp2w2 = 0.0; + + // initialize delay lines for pink filter network + st2095->lp1 = 0.0; + st2095->lp2 = 0.0; + st2095->lp3 = 0.0; + st2095->lp4 = 0.0; + st2095->lp5 = 0.0; + st2095->lp6 = 0.0; + + // cycle the generator for one complete time series to populate filter-bank delay lines + for (int i=0; isamplesPerPeriod; i++) + generate_st2095_noise_sample(st2095); + st2095->accum = 0.0; +} + +float generate_st2095_noise_sample( st2095_noise_t *st2095 ) { + float white, w, pink; + + // Generate a pseudorandom integer in the range 0 <= seed <= randMax. + //# Bitwise AND with randMax zeroes out any unwanted high order bits. + st2095->seed = (1664525 * st2095->seed + st2095->randStep) & st2095->randMax; + // Scale to a real number in the range -1.0 <= white <= 1.0 + white = (float)st2095->seed * st2095->scaleFactor - 1.0; + + // Run pink filter; a parallel network of first-order LP filters, scaled to + // produce an output signal with target RMS = -21.5 dB FS (-18.5 dB AES FS) + // when bandpass filter cutoff frequencies are 10 Hz and 22.4 kHz. + st2095->lp1 = 0.9994551 * st2095->lp1 + 0.00198166688621989 * white; + st2095->lp2 = 0.9969859 * st2095->lp2 + 0.00263702334184061 * white; + st2095->lp3 = 0.9844470 * st2095->lp3 + 0.00643213710202331 * white; + st2095->lp4 = 0.9161757 * st2095->lp4 + 0.01438952538362820 * white; + st2095->lp5 = 0.6563399 * st2095->lp5 + 0.02698408541064610 * white; + pink = st2095->lp1 + st2095->lp2 + st2095->lp3 + + st2095->lp4 + st2095->lp5 + st2095->lp6 + white * 0.0342675832159306; + st2095->lp6 = white * 0.0088766118009356; + + // Run bandpass filter; a series network of 4 biquad filters + // Biquad filters implemented in Direct Form II + w = pink - st2095->hp1_a1 * st2095->hp1w1 - st2095->hp1_a2 * st2095->hp1w2; + pink = st2095->hp1_b0 * w + st2095->hp1_b1 * st2095->hp1w1 + st2095->hp1_b2 * st2095->hp1w2; + st2095->hp1w2 = st2095->hp1w1; + st2095->hp1w1 = w; + + w = pink - st2095->hp2_a1 * st2095->hp2w1 - st2095->hp2_a2 * st2095->hp2w2; + pink = st2095->hp2_b0 * w + st2095->hp2_b1 * st2095->hp2w1 + st2095->hp2_b2 * st2095->hp2w2; + st2095->hp2w2 = st2095->hp2w1; + st2095->hp2w1 = w; + + w = pink - st2095->lp1_a1 * st2095->lp1w1 - st2095->lp1_a2 * st2095->lp1w2; + pink = st2095->lp1_b0 * w + st2095->lp1_b1 * st2095->lp1w1 + st2095->lp1_b2 * st2095->lp1w2; + st2095->lp1w2 = st2095->lp1w1; + st2095->lp1w1 = w; + + w = pink - st2095->lp2_a1 * st2095->lp2w1 - st2095->lp2_a2 * st2095->lp2w2; + pink = st2095->lp2_b0 * w + st2095->lp2_b1 * st2095->lp2w1 + st2095->lp2_b2 * st2095->lp2w2; + st2095->lp2w2 = st2095->lp2w1; + st2095->lp2w1 = w; + + // Limit peaks to +/-MaxAmp + if (pink > st2095->maxAmp) + pink = st2095->maxAmp; + else if (pink < -st2095->maxAmp) + pink = -st2095->maxAmp; + + // accumulate squared amplitude for RMS computation + st2095->accum += (pink * pink); + return(pink); +} diff --git a/speaker-test/st2095.h b/speaker-test/st2095.h new file mode 100644 index 0000000..32083da --- /dev/null +++ b/speaker-test/st2095.h @@ -0,0 +1,41 @@ +#define ST2095_MAX_PEAK -9.5 // dB +#define ST2095_HPFC 10.0 // Highpass filter cutoff in Hz +#define ST2095_LPFC 22400.0 // Lowpass filter cutoff in Hz + +typedef struct +{ + float maxAmp; + int samplesPerPeriod; + int randStep; + int randMax; + int seed; + float scaleFactor; + float w0t; + float k; + float k2; + float LpFc; + // biquad coefficients + float hp1_a1, hp1_a2; + float hp1_b0, hp1_b1, hp1_b2; + float hp2_a1, hp2_a2; + float hp2_b0, hp2_b1, hp2_b2; + float lp1_a1, lp1_a2; + float lp1_b0, lp1_b1, lp1_b2; + float lp2_a1, lp2_a2; + float lp2_b0, lp2_b1, lp2_b2; + // delay-line variables for bandpass filter + float hp1w1, hp1w2; + float hp2w1, hp2w2; + float lp1w1, lp1w2; + float lp2w1, lp2w2; + // delay-line variables for pink filter network + float lp1, lp2, lp3, lp4, lp5, lp6; + // statistics accumulator + float accum; +} st2095_noise_t; + +void initialize_st2095_noise( st2095_noise_t *st2095, int sample_rate ); +float generate_st2095_noise_sample( st2095_noise_t *st2095 ); + +void reset_st2095_noise_measurement( st2095_noise_t *st2095 ); +float compute_st2095_noise_measurement( st2095_noise_t *st2095, int period );