1/*
2  st2095.c
3
4  Generate Bandlimited Pink Noise (-18.5dB AES FS)
5  Using the SMPTE ST 2095:1-2015 standard
6
7  Based on pseudo-code from the above SMPTE standard, which bore the credit
8  "Revised 2015-01-04 by Calvert Dayton"
9
10  Copyleft 2023 Rick Sayre - No rights reserved.
11*/
12
13#include "aconfig.h"
14#include <stdio.h>
15#include <math.h>
16#include "st2095.h"
17
18/************************************************************/
19
20
21void reset_st2095_noise_measurement( st2095_noise_t *st2095 ) {
22    st2095->accum = 0.;
23}
24
25float compute_st2095_noise_measurement( st2095_noise_t *st2095, int period ) {
26    return(10. * log10f(st2095->accum / (float)period) + 3.01);
27}
28
29void initialize_st2095_noise( st2095_noise_t *st2095, int sample_rate) {
30    // Periodicity in samples must be a power of two, <= 2^31
31    // Typical values are 524288, 1048576, 2097152 or 4194304
32    if (sample_rate > 48000) {
33	// Special case LCG step for 1024K samples @ 88.2K or 96k
34	st2095->samplesPerPeriod = 1048576;
35	st2095->randStep = 163841;
36    } else {
37	st2095->samplesPerPeriod = 524288;
38	st2095->randStep = 52737;
39    }
40
41    // set up LCG PRNG
42    st2095->randMax = st2095->samplesPerPeriod - 1;
43    st2095->seed = 0;
44    st2095->scaleFactor = 2.0 / (float)st2095->randMax;
45
46    st2095->maxAmp = powf(10.0, ST2095_MAX_PEAK / 20.0);
47
48    // Calculate omegaT for matched Z transform highpass filters
49    st2095->w0t = 2.0 * M_PI * ST2095_HPFC / (float)sample_rate;
50
51    //  Limit LpFc <= Nyquist (actually lower, based on 48 vs 22.4 KHz spec cutoff)
52    //          The spec says the filter begins at 22.4KHz, if we ask for a Nyquist-impossible
53    //          sampling rate, compute something with the same relationship
54    st2095->LpFc = ST2095_LPFC;
55    float rateratio = 48000. / ST2095_LPFC;
56    if (st2095->LpFc > sample_rate/rateratio)
57	st2095->LpFc = sample_rate/rateratio;
58
59    // Calculate k and k^2 for bilinear transform lowpass filters
60    st2095->k = tanf(( 2.0 * M_PI * st2095->LpFc / (float)sample_rate ) / 2.0);
61    st2095->k2 = st2095->k * st2095->k;
62
63    // Calculate biquad coefficients for bandpass filter components
64    st2095->hp1_a1 = -2.0 * expf(-0.3826835 * st2095->w0t) * cosf(0.9238795 * st2095->w0t);
65    st2095->hp1_a2 = expf(2.0 * -0.3826835 * st2095->w0t);
66    st2095->hp1_b0 = (1.0 - st2095->hp1_a1 + st2095->hp1_a2) / 4.0;
67    st2095->hp1_b1 = -2.0 * st2095->hp1_b0;
68    st2095->hp1_b2 = st2095->hp1_b0;
69
70    st2095->hp2_a1 = -2.0 * expf(-0.9238795 * st2095->w0t) * cosf(0.3826835 * st2095->w0t);
71    st2095->hp2_a2 = expf(2.0 * -0.9238795 * st2095->w0t);
72    st2095->hp2_b0 = (1.0 - st2095->hp2_a1 + st2095->hp2_a2) / 4.0;
73    st2095->hp2_b1 = -2.0 * st2095->hp2_b0;
74    st2095->hp2_b2 = st2095->hp2_b0;
75
76    st2095->lp1_a1 = (2.0 * (st2095->k2 - 1.0)) / (st2095->k2 + (st2095->k / 1.306563) + 1.0);
77    st2095->lp1_a2 = (st2095->k2 - (st2095->k / 1.306563) + 1.0) / (st2095->k2 + (st2095->k / 1.306563) + 1.0);
78    st2095->lp1_b0 = st2095->k2 / (st2095->k2 + (st2095->k / 1.306563) + 1.0);
79    st2095->lp1_b1 = 2.0 * st2095->lp1_b0;
80    st2095->lp1_b2 = st2095->lp1_b0;
81
82    st2095->lp2_a1 = (2.0 * (st2095->k2 - 1.0)) / (st2095->k2 + (st2095->k / 0.541196) + 1.0);
83    st2095->lp2_a2 = (st2095->k2 - (st2095->k / 0.541196) + 1.0) / (st2095->k2 + (st2095->k / 0.541196) + 1.0);
84    st2095->lp2_b0 = st2095->k2 / (st2095->k2 + (st2095->k / 0.541196) + 1.0);
85    st2095->lp2_b1 = 2.0 * st2095->lp2_b0;
86    st2095->lp2_b2 = st2095->lp2_b0;
87
88    // initialize delay lines for bandpass filter
89    st2095->hp1w1 = 0.0;
90    st2095->hp1w2 = 0.0;
91    st2095->hp2w1 = 0.0;
92    st2095->hp2w2 = 0.0;
93    st2095->lp1w1 = 0.0;
94    st2095->lp1w2 = 0.0;
95    st2095->lp2w1 = 0.0;
96    st2095->lp2w2 = 0.0;
97
98    // initialize delay lines for pink filter network
99    st2095->lp1 = 0.0;
100    st2095->lp2 = 0.0;
101    st2095->lp3 = 0.0;
102    st2095->lp4 = 0.0;
103    st2095->lp5 = 0.0;
104    st2095->lp6 = 0.0;
105
106    // cycle the generator for one complete time series to populate filter-bank delay lines
107    for (int i=0; i<st2095->samplesPerPeriod; i++)
108	generate_st2095_noise_sample(st2095);
109    st2095->accum = 0.0;
110}
111
112float generate_st2095_noise_sample( st2095_noise_t *st2095 ) {
113    float white, w, pink;
114
115    // Generate a pseudorandom integer in the range 0 <= seed <= randMax.
116    //# Bitwise AND with randMax zeroes out any unwanted high order bits.
117    st2095->seed = (1664525 * st2095->seed + st2095->randStep) & st2095->randMax;
118    // Scale to a real number in the range -1.0 <= white <= 1.0
119    white = (float)st2095->seed * st2095->scaleFactor - 1.0;
120
121    // Run pink filter; a parallel network of first-order LP filters, scaled to
122    // produce an output signal with target RMS = -21.5 dB FS (-18.5 dB AES FS)
123    // when bandpass filter cutoff frequencies are 10 Hz and 22.4 kHz.
124    st2095->lp1 = 0.9994551 * st2095->lp1 + 0.00198166688621989 * white;
125    st2095->lp2 = 0.9969859 * st2095->lp2 + 0.00263702334184061 * white;
126    st2095->lp3 = 0.9844470 * st2095->lp3 + 0.00643213710202331 * white;
127    st2095->lp4 = 0.9161757 * st2095->lp4 + 0.01438952538362820 * white;
128    st2095->lp5 = 0.6563399 * st2095->lp5 + 0.02698408541064610 * white;
129    pink = st2095->lp1 + st2095->lp2 + st2095->lp3 +
130	st2095->lp4 + st2095->lp5 + st2095->lp6 + white * 0.0342675832159306;
131    st2095->lp6 = white * 0.0088766118009356;
132
133    // Run bandpass filter; a series network of 4 biquad filters
134    // Biquad filters implemented in Direct Form II
135    w = pink - st2095->hp1_a1 * st2095->hp1w1 - st2095->hp1_a2 * st2095->hp1w2;
136    pink = st2095->hp1_b0 * w + st2095->hp1_b1 * st2095->hp1w1 + st2095->hp1_b2 * st2095->hp1w2;
137    st2095->hp1w2 = st2095->hp1w1;
138    st2095->hp1w1 = w;
139
140    w = pink - st2095->hp2_a1 * st2095->hp2w1 - st2095->hp2_a2 * st2095->hp2w2;
141    pink = st2095->hp2_b0 * w + st2095->hp2_b1 * st2095->hp2w1 + st2095->hp2_b2 * st2095->hp2w2;
142    st2095->hp2w2 = st2095->hp2w1;
143    st2095->hp2w1 = w;
144
145    w = pink - st2095->lp1_a1 * st2095->lp1w1 - st2095->lp1_a2 * st2095->lp1w2;
146    pink = st2095->lp1_b0 * w + st2095->lp1_b1 * st2095->lp1w1 + st2095->lp1_b2 * st2095->lp1w2;
147    st2095->lp1w2 = st2095->lp1w1;
148    st2095->lp1w1 = w;
149
150    w = pink - st2095->lp2_a1 * st2095->lp2w1 - st2095->lp2_a2 * st2095->lp2w2;
151    pink = st2095->lp2_b0 * w + st2095->lp2_b1 * st2095->lp2w1 + st2095->lp2_b2 * st2095->lp2w2;
152    st2095->lp2w2 = st2095->lp2w1;
153    st2095->lp2w1 = w;
154
155    // Limit peaks to +/-MaxAmp
156    if (pink > st2095->maxAmp)
157	pink = st2095->maxAmp;
158    else if (pink < -st2095->maxAmp)
159	pink = -st2095->maxAmp;
160
161    // accumulate squared amplitude for RMS computation
162    st2095->accum += (pink * pink);
163    return(pink);
164}
165