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 <whorfin@gmail.com>
Signed-off-by: Jaroslav Kysela <perex@perex.cz>
This commit is contained in:
Rick Sayre 2023-12-30 19:11:48 -08:00 committed by Jaroslav Kysela
parent 6da3737bd4
commit e523020eb8
5 changed files with 272 additions and 19 deletions

View file

@ -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

View file

@ -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.

View file

@ -25,8 +25,13 @@
* Some cleanup from Daniel Caujolle-Bert <segfault@club-internet.fr>
* 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 <sys/time.h>
#include <math.h>
#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;

164
speaker-test/st2095.c Normal file
View file

@ -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 <stdio.h>
#include <math.h>
#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; i<st2095->samplesPerPeriod; 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);
}

41
speaker-test/st2095.h Normal file
View file

@ -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 );