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.
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
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");
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;
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
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.
The loudness is 42.5 sones, with a peak at 6 (Bark scale). Convert the peak loudness value to Hz using
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.
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
The perceived sharpness of a sound might contribute significantly to how annoying it is to people. Estimate the perceived sharpness level using
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.
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).
acousticFluctuation to compute the perceived fluctuation strength. Let the algorithm automatically detect the most audible fluctuation frequency (fMod).
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
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:
A quantitative description was developed using the results of psychoacoustic experiments:
percentile loudness in sone (level that is exceeded only 5% of the time)
for , where is the sharpness in acum
, where is the fluctuation strength in vacil and 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 and to zero.
Percentile loudness, , is the second value returned by the third output of
"TimeVarying" is set to
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
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
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
[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.
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
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