Main Content

Effect of Hearing Protection on Perceived Noise Levels

In this example, you measure the noise of a helicopter and use psychoacoustic metrics to model its perceived loudness, sharpness, fluctuation strength, and relative annoyance level. You then model the effect of hearing protection to estimate the perceived improvement.

Recording Level Calibration

Measuring sound pressure levels (SPL) or perceptual metrics such as loudness requires you to first calibrate your microphone input level. You can use calibrateMicrophone to match your recording level to the SPL reading of a calibrated meter. All you need is a calibrated SPL meter and a 1 kHz tone being played back from your computer or even from an app on your phone. The 1 kHz calibration tone does not have to be at a known calibrated level, but it should be loud enough to dominate any background noise. For an example of how to calibrate by playing the tone with MATLAB, see the documentation for calibrateMicrophone.

Simulate the tone recording and include some background noise. Assume an SPL meter gave a reading of 90.1 dB (C-weighted).

FS = 48e3;
t = (1:2*FS)/FS;
s = rng('default');
xtone = 0.046*sin(2*pi*t*1000).' + 1e-2*pinknoise(2*FS);
rng(s)

splMeterReading = 90.1; 

To compute the calibration level of a recording chain, call calibrateMicrophone and specify the test tone, the sample rate, the SPL reading, and the frequency weighting of the SPL meter. Even though the signal at 1 kHz is not affected by the frequency weighting, the background noise is, so you get a more precise calibration level by matching the meter's frequency weighting.

calib = calibrateMicrophone(xtone,FS,splMeterReading,"FrequencyWeighting","C-weighting");

Sound Pressure Levels (SPL)

Once you have a calibration factor for a recording chain, you can make acoustic measurements. When using a physical meter, you might be limited to whatever settings were selected at the time of the measurement. Here, using your recording and the splMeter object, you can select any settings. This makes it easy to experiment with different time and frequency weighting options.

Load a helicopter recording that contains spatial information. You only need to keep the first channel, which corresponds to the omnidirectional signal (similar to a mono recording with a normal omnidirectional microphone). Create the corresponding time vector.

[x,FS] = audioread("Heli_16ch_ACN_SN3D.wav");
x = x(:,1);
t = (1:size(x,1))/FS;

Create an splMeter object and select C-weighting, fast time weighting, and a 0.2 second interval for the peak SPL measurement.

spl = splMeter("CalibrationFactor",calib, ...
               "FrequencyWeighting","C-weighting", ...
               "TimeWeighting","Fast", ...
               "TimeInterval",0.2, ...
               "SampleRate",FS);

Plot SPL and peak SPL.

[LCF,~,LCpeak] = spl(x);
plot(t,LCpeak,t,LCF)
legend('LCpeak','LCF','Location',"southeast")
title('SPL Measurement of Helicopter Noise')
xlabel('Time (seconds)')
ylabel('SPL (dB)')
ylim([68 100])
grid on

Psychoacoustic Metrics

Loudness level

Monitoring SPL is important for hearing safety, but acousticLoudness allows you to estimate loudness levels as they are actually perceived by people with normal hearing (that is, no large hearing impairments). You can also see what frequency bands contribute the most to the sensation of loudness.

Using the same calibration level as before, and assuming a free-field recording (the default), plot stationary loudness.

acousticLoudness(x,FS,calib)

The loudness is 42.5 sones, with a peak at 6 (Bark scale). Convert the peak loudness value to Hz using bark2hz

fprintf("Peak loudness frequency: %d Hz\n",round(bark2hz(6),-1));
Peak loudness frequency: 630 Hz

By default, the levels are in sones. It might be easier to understand this measurement if it was compared to an SPL (dB) reading. The phons unit of loudness is matched to a reference frequency of 1 kHz. For example, any signal with a loudness of 60 phons is perceived to be as loud as a 1 kHz tone that registers at 60 dB on an SPL meter (1 kHz is unaffected by the frequency weighting). Consequently, converting 42.5 sones to phons using sone2phon lets you know that the helicopter is percieved as loud as a 94 dB 1 kHz tone.

fprintf("Equivalent 1 kHz SPL: %d phons\n", round(sone2phon(42.5)));
Equivalent 1 kHz SPL: 94 phons

Make your own plot with units in phons and frequency in Hz on a log scale.

[sone,spec] = acousticLoudness(x,FS,calib);
barks = 0.1:0.1:24; % Bark scale for ISO 532-1 loudness
hz = bark2hz(barks);
specPhon = sone2phon(spec);
semilogx(hz,specPhon)
title('Specific Loudness')
subtitle(sprintf('Loudness = %.1f phons',sone2phon(sone)))
xlabel('Frequency (Hz)')
ylabel('Loudness (phons/Bark)')
xlim(hz([1,end]))
grid on

You can also plot time-varying loudness and specific loudness to analyze the sound of the moving helicopter. In this case, the noise is stationary, but you can observe the distictive "impulsive" nature of a helicopter.

acousticLoudness(x,FS,calib,'TimeVarying',true)

Plot specific loudness with the frequency in Hz (log-scale).

[tvsone,tvspec,perc] = acousticLoudness(x,FS,calib,'TimeVarying',true);
spectime = 0:2e-3:2e-3*(size(tvspec,1)-1); % time axis for loudness output
clf; % do not reuse the previous subplot
surf(spectime,hz,sone2phon(tvspec).','EdgeColor','interp');
set(gca,'View',[0 90],'YScale','log','YLim',hz([1,end]));
title('Specific Loudness')
zlabel('Specific Loudness (phons/Bark)')
ylabel('Frequency (Hz)')
xlabel('Time (seconds)')
colorbar

Sharpness Level

The perceived sharpness of a sound might contribute significantly to how annoying it is to people. Estimate the perceived sharpness level using acousticSharpness.

sharp = acousticSharpness(spec)
sharp = 1.4507

Pink noise has a sharpness of 2 acums, so you now know that the helicopter noise is biased towards low frequencies.

Plot time-varying sharpness.

acousticSharpness(x,FS,calib,'TimeVarying',true);

Fluctuation Strength

In the case of helicopter noise, the low-frequency modulation also contributes to how annoying the noise might be perceived.

First, look at what modulation frequencies are found in the signal.

N = 2^nextpow2(size(x,1));
xa = abs(x); % Use the rectified signal
pspectrum(xa-mean(xa),FS,'FrequencyLimits',[0 80],'FrequencyResolution',1)
title('Modulation Frequencies')

Observe that the modulation frequency peaks at 18.8 Hz. Below 30 Hz, modulation is perceived as fluctuation (as opposed to roughness).

Use acousticFluctuation to compute the perceived fluctuation strength. Let the algorithm automatically detect the most audible fluctuation frequency (fMod).

acousticFluctuation(x,FS,calib)

You might want to use Hz instead of Bark. To save computations, reuse the time-varying specific loudness that you computed previously.

[vacil,spec,fMod] = acousticFluctuation(tvspec);
clf; % do not reuse previous subplot
flucHz = bark2hz(0.5:0.5:23.5);
surf(spectime,flucHz,spec.','EdgeColor','interp');
set(gca,'View',[0 90],'YScale','log','YLim',flucHz([1,end]));
title('Specific Flucutation Strength')
zlabel('Specific Flucutation Strength (vacils/Bark)')
ylabel('Frequency (Hz)')
xlabel('Time (seconds)')
colorbar

Sound Quality

For the evaluation of sound quality, the previous metrics can be combined to produce a psychoacoustic annoyance factor (Zwicker and Fastl). The relation is as follows:

PAN(1+[g1(S)]2+[g2(F,R)]2)

A quantitative description was developed using the results of psychoacoustic experiments:

PA=N5(1+wS2+wFR2)

with:

  • N5 percentile loudness in sone (level that is exceeded only 5% of the time)

  • wS=(S-1.75)(0.25log10(N5+10)) for S>1.75, where S is the sharpness in acum

  • wFR=2.18(N5)0.4(0.4F+0.6R), where F is the fluctuation strength in vacil and R is the roughness in asper

In this example, sharpness was less than 1.75, so it is not a contributing factor. The noise is mainly modulated by a frequency lower than 30 Hz, so consider roughness negligible for this example. Therefore, you can set R and ws to zero.

Percentile loudness, N5, is the second value returned by the third output of acousticLoudness when "TimeVarying" is set to true.

N5 = perc(2);

Compute an average fluctuation strength ignoring the first second of the signal.

f = mean(vacil(501:end,1));

Compute the psychoacoustic annoyance factor.

pa = N5 * (1 + abs(2.18/(N5^.4)*(.4*f)))
pa = 46.7171

Effect of Ear Protection

Measure the impact of hearing protection on the measured SPL and the perceived noise. The attenuation characteristics of a 24-dB-rated hearing protector are:

Frequency (Hz) 125 250 500 1000 2000 3000 4000 6000 8000

Attenuation (dB) 14.4 21.6 27.0 31.9 36.1 39.5 41.9 39.8 37.0

Interpolate the hearing protection data.

att = [  125   250   500  1000  2000  3000  4000  6000  8000 ; ...
        14.4  21.6  27.0  31.9  36.1  39.5  41.9  39.8  37.0 ].';

att = [ 20,0; att ]; % Assume zero attenuation at 20 Hz

cf = getCenterFrequencies(splMeter("Bandwidth","1/3 octave"));

p = interp1(att(:,1),att(:,2),cf,"pchip");

Display the interpolated attenuation.

semilogx(cf,-p,att(:,1),-att(:,2),'rx')
title('Attenuation from Hearing Protection')
ylabel('Relative SPL (dB)')
xlabel('Frequency (Hz)')
xlim(cf([1,end]))
grid on

Simulation Using Graphic EQ Filter Bank

Design a graphicEQ object to simulate the effect of the hearing protection.

geq = graphicEQ("SampleRate",FS,"Bandwidth","1/3 octave","EQOrder",14,"Gains",-p);

Display the frequency response of the graphicEQ object.

[B,A] = coeffs(geq);
sos = [B;[ones(1,size(A,2));A]].';
[H,w] = freqz(sos,2^16,FS);
semilogx(w,db(abs(H)))
title('Frequency Response of Hearing Protection Simulation Filter')
ylabel('Relative SPL (dB)')
xlabel('Frequency (Hz)')
xlim(cf([1,end]))
grid on

Filter the helicoptor recording using the graphic EQ to simulate the hearing protection.

xhp = geq(x);

Compare the SPL with and without hearing protection. Reuse the same SPL meter settings, but use the filtered recording. Set the Y limits as before to make visual comparison easier with the earlier plot.

reset(spl)
[LCFnew,~,LCpeaknew] = spl(xhp);
plot(t,LCpeak,t,LCF,t,LCpeaknew,t,LCFnew)
legend('LCpeak (original)', 'LCF (original)', ...
       'LCpeak (with ear protection)', ...
       'LCF (with ear protection)', ...
       'Location','southeast')
title('SPL Measurement of Helicopter Noise')
xlabel('Time (seconds)')
ylabel('SPL (dB)')
ylim([68 100])
grid on

The attenuation might be less than expected by looking at the specification, but as you saw with the sharpness measurement, the heliocopter noise is concentrated in the lower frequencies where any hearing protection is less efficient.

Look at the perceived loudness to get a better idea of how loud the helicopter now sounds.

acousticLoudness(xhp,FS,calib)

Loudness went from 42.5 to 6.2 sones. However, it might be easier to look at loudness in phons instead of sones, because phons are equivalent to the SPL of a 1 kHz tone. Convert the sone units to phons using sone2phon.

fprintf("Loudness without hearing protection:\t%.1f phons\n",sone2phon(42.5));
Loudness without hearing protection:	94.1 phons
fprintf("Loudness with hearing protection:\t%.1f phons\n",sone2phon(6.2));
Loudness with hearing protection:	66.3 phons
fprintf("Perceived noise reduction:\t\t%.1f phons (dB SPL at 1 kHz)\n",sone2phon(42.5)-sone2phon(6.2));
Perceived noise reduction:		27.8 phons (dB SPL at 1 kHz)

This hearing protection is rated as 24 dB, but the helicopter noise is perceived as approximately 28 dB quieter when wearing them.

Calculate the reduction in the psychoacoustic annoyance factor. Start by computing the mean sharpness.

[~,spechp,perchp] = acousticLoudness(xhp,FS,calib,"TimeVarying",true);
shp = acousticSharpness(spechp,'TimeVarying',true);
new_sharp = mean(shp(501:end))
new_sharp = 0.7258

Sharpness is below the threshold of 1.75, so it is ignored. Now, compute the mean fluctuation strength.

vacilhp = acousticFluctuation(spechp);
fhp = mean(vacilhp(501:end,1));

Compute the new psychoacoustic annoyance factor. It has reduced from approximately 47 to 7.

N5hp = perchp(2); % N5 with hearing protection
pahp = N5hp * (1 + abs(2.18/(N5hp^.4)*(.4*fhp)))
pahp = 7.0230