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 
reset_st2095_noise_measurement( st2095_noise_t *st2095 )21 void reset_st2095_noise_measurement( st2095_noise_t *st2095 ) {
22     st2095->accum = 0.;
23 }
24 
compute_st2095_noise_measurement( st2095_noise_t *st2095, int period )25 float compute_st2095_noise_measurement( st2095_noise_t *st2095, int period ) {
26     return(10. * log10f(st2095->accum / (float)period) + 3.01);
27 }
28 
initialize_st2095_noise( st2095_noise_t *st2095, int sample_rate)29 void 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 
generate_st2095_noise_sample( st2095_noise_t *st2095 )112 float 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