Acoustics engineer: MQL recipes

Recipes

Adding noise sources and exposure time
Doppler shift of a passing siren
Spectrum and THD of a tone
Reverberation time of a classroom
Fan noise PSD: periodogram against Welch
Tuning a bass-reflex port as a Helmholtz resonator
Room modes and degenerate frequencies of a listening room
Sabine against Eyring in a dead recording booth
Noise barrier: energy-average Leq and a significance test
Absorption of glass wool with the Delany-Bazley model

executed on MentDB

Adding noise sources and exposure time

Three machines measured alone at 1 m give 85, 88 and 90 dB. Add them on an energy basis, find the distance where the level falls to 80 dB (free field, -20 log r), and the daily exposure allowed at 1 m (85 dB for 8 h, 3 dB exchange rate).

Pipeline

01 · sumL = 10 log10(sum 10^(Li/10))math eval→02 · distanceSolve L - 20 log10(r) = 80num solve→03 · exposureT = 8 h / 2^((L-85)/3)math calc

Script

admin

#1. Three machines measured alone at 1 m: 85, 88 and 90 dB;
-> "[L]" (math eval "10*log10(10^(a/10)+10^(b/10)+10^(c/10))" "{\"a\": 85, \"b\": 88, \"c\": 90}");
#2. Free-field point source: -20 log10 r2/r1, distance where the level falls to 80 dB;
-> "[r]" (num solve (concat [L] "-20*log10(r)-80") "r" 1 100);
#3. Daily exposure allowed at 1 m, 85 dB criterion and 3 dB exchange rate;
-> "[T]" (math calc "8/2^(([L]-85)/3)");
concat "total = " (math round [L] 1) " dB, 80 dB at " (math round [r] 2) " m; allowed exposure at 1 m = " (math round (math calc "[T]*60") 0) " min per day"

Result

mentdb

total = 92.9 dB, 80 dB at 4.41 m; allowed exposure at 1 m = 77 min per day

Go further

Add the air absorption term with num solve for long distances
Compute a daily noise dose over several work stations with num sum

executed on MentDB

Doppler shift of a passing siren

An ambulance drives at 90 km/h with a 700 Hz siren in 25 degC air. Compute the speed of sound, the frequency heard while it approaches and while it moves away, and the pitch drop in semitones.

Pipeline

01 · soundc = 331.3 sqrt(1 + T/273.15)math calc→02 · speedSpeed in m/sunit convert→03 · dopplerf c/(c - v) and f c/(c + v)math calc→04 · pitch12 log2(fa/fr)math round

Script

admin

#1. Speed of sound at 25 degC;
-> "[c]" (math calc "331.3*sqrt(1+25/273.15)");
#2. Ambulance at 90 km/h, siren at 700 Hz;
-> "[v]" (unit convert 90 "km/h" "m/s");
-> "[fa]" (math calc "700*[c]/([c]-[v])");
-> "[fr]" (math calc "700*[c]/([c]+[v])");
#3. Pitch drop heard when it passes, in semitones 12 log2 fa/fr;
concat "c = " (math round [c] 1) " m/s; approaching " (math round [fa] 1) " Hz, receding " (math round [fr] 1) " Hz: a drop of " (math round (math calc "12*log_base([fa]/[fr], 2)") 2) " semitones"

Result

mentdb

c = 346.1 m/s; approaching 754.5 Hz, receding 652.8 Hz: a drop of 2.51 semitones

Go further

Find the speed that gives a drop of one full tone with num solve
Add a moving listener: f (c + vo)/(c - vs)

executed on MentDB

Spectrum and THD of a tone

Sample a 200 Hz tone with a 2nd and a 3rd harmonic at 4 kHz (200 samples), compute its amplitude spectrum, detect the peaks and the total harmonic distortion THD = sqrt(A2^2 + A3^2)/A1.

Pipeline

01 · sampleSamples of the signalnum table→02 · columnKeep the valuesmatrix col→03 · spectrumAmplitude spectrumdsp spectrum→04 · peaksHarmonic peaksdsp peaks

Script

admin

#1. 200 Hz tone with 2nd and 3rd harmonics, sampled at 4 kHz, 200 samples;
-> "[x]" (matrix col (num table "sin(2*pi*200*t)+0.3*sin(2*pi*400*t)+0.1*sin(2*pi*600*t)" "t" 0 0.04975 0.00025) 2);
#2. Amplitude spectrum, resolution fs/n = 20 Hz;
json load "rx_sp" (dsp spectrum [x] 4000);
#3. Peaks above 0.05;
json load "rx_pk" (dsp peaks (json select "rx_sp" /amplitude) "{\"height\": 0.05}");
-> "[A1]" (json select "rx_pk" /values[0]); -> "[A2]" (json select "rx_pk" /values[1]); -> "[A3]" (json select "rx_pk" /values[2]);
#4. Total harmonic distortion;
concat "peaks at " (math calc (concat "20*" (json select "rx_pk" /indices[0]))) ", " (math calc (concat "20*" (json select "rx_pk" /indices[1]))) ", " (math calc (concat "20*" (json select "rx_pk" /indices[2]))) " Hz, amplitudes " (math round [A1] 4) ", " (math round [A2] 4) ", " (math round [A3] 4) "; THD = " (math round (math calc "sqrt([A2]^2+[A3]^2)/[A1]*100") 2) " %"

Result

mentdb

peaks at 200, 400, 600 Hz, amplitudes 1, 0.3, 0.1; THD = 31.62 %

Go further

Use 201 samples to see the leakage, then reduce it with dsp window "hann"
Draw the spectrum with plot spectrum

executed on MentDB

Reverberation time of a classroom

A 10 x 8 x 3.2 m classroom has a tiled floor, an acoustic ceiling, painted walls and windows. Compute the absorption area at 500 Hz, the Sabine and Eyring RT60, and the area of absorbing panels needed on the walls to reach 0.5 s.

Pipeline

01 · absorptionA = sum S_i alpha_ivector dot→02 · sabineRT60 = 0.161 V / Amath calc→03 · eyringEyring with the total surfacestat sum→04 · treatmentPanels for the targetmath round

Script

admin

#1. Classroom 10 x 8 x 3.2 m: floor, ceiling tiles, walls, windows in m2;
-> "[S]" "[80, 80, 103.2, 12]";
-> "[alpha]" "[0.05, 0.70, 0.03, 0.18]";
#2. Equivalent absorption area at 500 Hz;
-> "[A]" (vector dot [S] [alpha]);
-> "[V]" (math calc "10*8*3.2");
#3. Sabine and Eyring reverberation times;
-> "[tsab]" (math calc "0.161*[V]/[A]");
-> "[am]" (math calc (concat "[A]/" (stat sum [S])));
-> "[teyr]" (math calc (concat "0.161*[V]/(-" (stat sum [S]) "*ln(1-[am]))"));
#4. Absorbing panels alpha 0.90 on the walls to reach 0.5 s;
concat "A = " (math round [A] 2) " m2, RT60 Sabine = " (math round [tsab] 2) " s, Eyring = " (math round [teyr] 2) " s; for 0.5 s add " (math round (math calc "(0.161*[V]/0.5-[A])/(0.90-0.03)") 1) " m2 of panels"

Result

mentdb

A = 65.26 m2, RT60 Sabine = 0.63 s, Eyring = 0.55 s; for 0.5 s add 19.7 m2 of panels

Go further

Repeat for each octave band with a matrix of alphas and matrix mul
Compute the critical distance sqrt(A/(16 pi))

executed on MentDB

Fan noise PSD: periodogram against Welch

A fan noise record (blade passing tone at 160 Hz, harmonic at 480 Hz, broadband noise) is analysed like in a vibration analyser: Hann window factors, raw periodogram and Welch averaged PSD, drawn together on a log scale.

Pipeline

01 · windowHann window, CG and ENBWdsp window→02 · rawPeriodogramdsp periodogram→03 · averageWelch PSDdsp welch→04 · drawOverlay of both PSDplot chart

Script

admin

#1. Fan noise sampled at 2 kHz, 2048 values: blade passing tone 160 Hz, harmonic 480 Hz, broadband noise;
-> "[fs]" 2000;
-> "[x]" (vector add (matrix col (num table "sin(2*pi*160*t)+0.5*sin(2*pi*480*t)" "t" 0 1.0235 0.0005) 2) (stat dist_random "normal" "[0, 1]" 2048));
#2. Hann window of 256 points: coherent gain and equivalent noise bandwidth;
-> "[w]" (dsp window "hann" 256);
-> "[cg]" (math round (stat mean [w]) 3);
-> "[enbw]" (math round (math calc (concat "256*" (stat sum_sq [w]) "/" (stat sum [w]) "^2*[fs]/256")) 2);
#3. Raw periodogram and Welch average of 50 % overlapping Hann segments;
json load "rx_per" (dsp periodogram [x] [fs]);
json load "rx_wel" (dsp welch [x] [fs] 256);
#4. Overlay both power spectral densities, log scale;
plot chart (concat "{\"title\": \"Fan noise PSD - Welch " (json select "rx_wel" /segments) " segments, Hann CG " [cg] ", ENBW " [enbw] " Hz\", \"x_label\": \"Frequency (Hz)\", \"y_label\": \"PSD (Pa2/Hz)\", \"log_y\": true, \"series\": [{\"name\": \"periodogram\", \"x\": " (json select "rx_per" /f) ", \"y\": " (json select "rx_per" /psd) "}, {\"name\": \"Welch\", \"x\": " (json select "rx_wel" /f) ", \"y\": " (json select "rx_wel" /psd) "}]}")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Find the tones in the Welch PSD with dsp peaks and a prominence condition.
Compare the windows: dsp window "flattop" 256 gives the best amplitude accuracy.

executed on MentDB

Tuning a bass-reflex port as a Helmholtz resonator

A 40 L loudspeaker box has a round port 7 cm wide and 15 cm long. The script adds the end corrections, gives the Helmholtz tuning frequency and the port length that tunes the box to 38 Hz.

Pipeline

01 · convertBox volume L to m3unit convert→02 · constant2 pi from the constant familyconstant math_two_pi→03 · parseValue of the constantatom get→04 · computefb = c/2pi sqrt(S/V Le)math calc

Script

admin

#1. Bass-reflex box: 40 L, round port 7 cm wide and 15 cm long, air at 20 degC;
-> "[V]" (unit convert 40 "L" "m3");
-> "[c]" (math calc "331.3*sqrt(1+20/273.15)");
-> "[r]" 0.035; -> "[S]" (math calc "pi*[r]^2");
#2. Effective length with the end corrections, 0.85 r inside + 0.85 r outside;
-> "[Le]" (math calc "0.15+1.7*[r]");
#3. Helmholtz frequency f = c / 2 pi x sqrt S / V Le, 2 pi from the constant family;
-> "[twopi]" (string replace (atom get (constant math_two_pi) 1 " ") "(" "");
-> "[fb]" (math calc "[c]/[twopi]*sqrt([S]/([V]*[Le]))");
#4. Port length needed to tune the box to 38 Hz;
-> "[L38]" (math calc "[S]*[c]^2/(([twopi]*38)^2*[V])-1.7*[r]");
concat "V = " [V] " m3, c = " (math round [c] 1) " m/s, Le = " (math round (math calc "[Le]*100") 2) " cm: fb = " (math round [fb] 1) " Hz; for fb = 38 Hz cut the port to " (math round (math calc "[L38]*100") 1) " cm"

Result

mentdb

V = 0.04 m3, c = 343.2 m/s, Le = 20.95 cm: fb = 37 Hz; for fb = 38 Hz cut the port to 13.9 cm

Go further

Use num solve to find the port diameter that keeps the air speed below 17 m/s at 38 Hz
Compare with a closed box: resonance of the driver with the air spring of the box

executed on MentDB

Room modes and degenerate frequencies of a listening room

A 6 x 4 x 3 m room has bad proportions. The script lists the axial, tangential and oblique modes up to 120 Hz, sorts them and finds the frequencies where several modes pile up, then gives the Schroeder frequency.

Pipeline

01 · modesf = c/2 sqrt(sum (n/L)^2) with nested hypotmath hypot→02 · roundNearest hertzmath rint→03 · sortSorted mode listmath order_double→04 · degenerateMost repeated frequencystat mode→05 · schroederfs = 2000 sqrt(T/V)math calc

Script

admin

#1. Rectangular listening room 6.0 x 4.0 x 3.0 m, c = 343 m/s;
-> "[c]" 343; -> "[Lx]" 6.0; -> "[Ly]" 4.0; -> "[Lz]" 3.0;
-> "[f]" "";
#2. Modes f = c/2 x sqrt of nx/Lx^2 + ny/Ly^2 + nz/Lz^2, with two nested hypot;
for (-> "[nx]" 0) (<= [nx] 2) (++ "[nx]") {
	for (-> "[ny]" 0) (<= [ny] 2) (++ "[ny]") {
		for (-> "[nz]" 0) (<= [nz] 1) (++ "[nz]") {
			if (> (math calc "[nx]+[ny]+[nz]") 0) {
				-> "[fm]" (math calc (concat [c] "/2*" (math hypot (math hypot (math calc "[nx]/[Lx]") (math calc "[ny]/[Ly]")) (math calc "[nz]/[Lz]"))));
				-> "[f]" (concat [f] "," (math rint [fm]));
			};
		};
	};
};
#3. Sorted mode list and the most repeated frequency = degenerate modes;
-> "[sorted]" (math order_double (concat "[" (string substr [f] 1) "]"));
-> "[deg]" (stat mode [sorted]);
#4. Schroeder frequency for T = 0.4 s;
-> "[fs]" (math calc "2000*sqrt(0.4/([Lx]*[Ly]*[Lz]))");
concat "modes (Hz): " [sorted] "; degenerate: " [deg] " Hz; Schroeder frequency " (math round [fs] 0) " Hz"

Result

mentdb

modes (Hz): [29.0,43.0,52.0,57.0,57.0,64.0,71.0,71.0,77.0,81.0,86.0,90.0,92.0,103.0,103.0,107.0,118.0]; degenerate: [57,71,103] Hz; Schroeder frequency 149 Hz

Go further

Try the Bolt golden ratios 1 : 1.28 : 1.54 and check that stat mode finds fewer coincidences
Count the modes per third-octave band with stat histogram

executed on MentDB

Sabine against Eyring in a dead recording booth

In a small, very absorbing booth Sabine's formula overestimates the reverberation time. The script computes the area-weighted absorption, the Sabine and Eyring times and the critical distance.

Pipeline

01 · meanArea-weighted mean absorptionstat weighted_mean→02 · absorptionTotal absorption A = sum S alphavector dot→03 · eyring-S ln(1 - alpha) with log1pmath log1p→04 · distanceCritical distancemath calc

Script

admin

#1. Recording booth 4 x 3 x 2.5 m: floor, ceiling, walls with their absorption at 1 kHz;
-> "[V]" (math calc "4*3*2.5");
-> "[areas]" "[12, 12, 35]";
-> "[alpha]" "[0.10, 0.70, 0.55]";
#2. Area-weighted mean absorption and total absorption A;
-> "[S]" (stat sum [areas]);
-> "[am]" (stat weighted_mean [alpha] [areas]);
-> "[A]" (vector dot [areas] [alpha]);
#3. Sabine T = 0.161 V / A, Eyring T = 0.161 V / -S ln 1 - am, with log1p for accuracy;
-> "[Ts]" (math calc "0.161*[V]/[A]");
-> "[Te]" (math calc (concat "0.161*[V]/(-[S]*" (math log1p (math calc "-[am]")) ")"));
#4. Critical distance where direct and reverberant fields are equal, Q = 2;
-> "[rc]" (math calc "0.057*sqrt(2*[V]/[Te])");
concat "mean alpha = " (math round [am] 3) ", A = " (math round [A] 2) " m2 Sabine; T Sabine = " (math round [Ts] 3) " s, T Eyring = " (math round [Te] 3) " s, Sabine overestimates by " (math round (math calc "([Ts]/[Te]-1)*100") 1) " %; critical distance " (math round [rc] 2) " m"

Result

mentdb

mean alpha = 0.489, A = 28.85 m2 Sabine; T Sabine = 0.167 s, T Eyring = 0.122 s, Sabine overestimates by 37.3 %; critical distance 1.26 m

Go further

Add air absorption 4 m V at high frequencies
Solve the extra absorption needed for T = 0.25 s with num solve

executed on MentDB

Noise barrier: energy-average Leq and a significance test

Hourly LAeq values were measured at a facade before and after a noise barrier. The script computes the energy-average Leq (not the arithmetic mean), the insertion loss and a Mann-Whitney test on the two samples.

Pipeline

01 · energy10^(L/10) of each hourmath calc→02 · averageMean of the energiesmath avg→03 · testNon-parametric comparisonstat mann_whitney→04 · compareArithmetic meansstat mean

Script

admin

#1. Hourly LAeq at the facade, dB, before and after the noise barrier;
-> "[before]" "68.2,70.1,71.4,69.8,72.3,70.6,69.1,71.0";
-> "[after]" "63.5,64.8,66.2,63.9,65.7,64.1,66.0";
#2. Energy average Leq = 10 log10 of the mean of 10^L/10, not the arithmetic mean;
-> "[eb]" ""; -> "[ea]" "";
for (-> "[i]" 1) (<= [i] (atom size [before] ",")) (++ "[i]") { -> "[eb]" (concat [eb] "," (math calc (concat "10^(" (atom get [before] [i] ",") "/10)"))); };
for (-> "[i]" 1) (<= [i] (atom size [after] ",")) (++ "[i]") { -> "[ea]" (concat [ea] "," (math calc (concat "10^(" (atom get [after] [i] ",") "/10)"))); };
-> "[Lb]" (math calc (concat "10*log10(" (math avg (concat "[" (string substr [eb] 1) "]")) ")"));
-> "[La]" (math calc (concat "10*log10(" (math avg (concat "[" (string substr [ea] 1) "]")) ")"));
#3. Is the reduction significant? Non-parametric test, hourly levels are not normal;
-> "[mw]" (stat mann_whitney (concat "[" [before] "]") (concat "[" [after] "]"));
concat "Leq before " (math round [Lb] 1) " dB(A), after " (math round [La] 1) " dB(A): insertion loss " (math round (math calc "[Lb]-[La]") 1) " dB; arithmetic means " (math round (stat mean (concat "[" [before] "]")) 1) " and " (math round (stat mean (concat "[" [after] "]")) 1) "; Mann-Whitney " [mw]

Result

mentdb

Leq before 70.5 dB(A), after 65 dB(A): insertion loss 5.5 dB; arithmetic means 70.3 and 64.9; Mann-Whitney {"u": 56, "p_value": 3.108003108003108E-4, "method": "exact"}

Go further

Compute L10 and L90 with stat percentile
Build the day-evening-night level Lden with the 5 and 10 dB penalties

executed on MentDB

Absorption of glass wool with the Delany-Bazley model

A 50 mm layer of glass wool lies on a rigid wall. With the Delany-Bazley model, the script computes the complex characteristic impedance and wave number, the surface impedance -i Zc cot(kd) and the absorption coefficient at four octave bands.

Pipeline

01 · modelComplex Zc and kcomplex sub→02 · surfaceZs = -i Zc / tan(kd)complex tan→03 · reflectionR = (Zs - Z0)/(Zs + Z0)complex div→04 · absorptionalpha = 1 - |R|^2complex abs

Script

admin

#1. Glass wool on a rigid wall: flow resistivity 20000 Pa.s/m2, thickness 50 mm;
-> "[sig]" 20000; -> "[d]" 0.05; -> "[rho]" 1.204; -> "[c]" 343.2; -> "[Z0]" (math calc "[rho]*[c]");
-> "[out]" "";
for (-> "[i]" 1) (<= [i] 4) (++ "[i]") {
	-> "[f]" (atom get "250,500,1000,2000" [i] ",");
	#2. Delany-Bazley model: complex characteristic impedance Zc and wave number k;
	-> "[X]" (math calc "[rho]*[f]/[sig]");
	-> "[Zc]" (complex sub (math calc "[Z0]*(1+0.0571*[X]^-0.754)") (concat (math calc "[Z0]*0.087*[X]^-0.732") "i"));
	-> "[k]" (complex sub (math calc "2*pi*[f]/[c]*(1+0.0978*[X]^-0.700)") (concat (math calc "2*pi*[f]/[c]*0.189*[X]^-0.595") "i"));
	#3. Surface impedance of the backed layer Zs = -i Zc cot kd;
	-> "[Zs]" (complex div (complex mul "-i" [Zc]) (complex tan (complex mul [k] [d])));
	#4. Reflection factor and absorption coefficient at normal incidence;
	-> "[R]" (complex div (complex sub [Zs] [Z0]) (complex add [Zs] [Z0]));
	-> "[out]" (concat [out] [f] " Hz: alpha = " (math round (math calc (concat "1-" (complex abs [R]) "^2")) 3) "; ");
};
[out]

Result

mentdb

250 Hz: alpha = 0.222; 500 Hz: alpha = 0.619; 1000 Hz: alpha = 0.943; 2000 Hz: alpha = 0.958;

Go further

Add an air gap behind the layer: transfer the impedance through the gap with a second complex tan
Plot alpha against frequency with plot chart on log_x