# MATLAB code to evaluate audio loundness over time

can anyone tell me how to evaluate the loudness/intensity over time using MATLAB? By loudness/intensity over time I mean this:

In the images above, the "black" one on the top is a sound's amplitude over time, and the green line is the corresponding loudness/intensity over time. The green line underneath, which is what I want, is generated using Praat.

Can any one tell me, that given a sound input, how can I generate the corresponding loudness/intensity over time using MATLAB?

For example, given a input, I can generate the upper one using this MATLAB script:

[x,fs] = audioread('soundfile.wav');
t = (0:length(x)-1)/fs;
subplot(211);
plot(t,x);


How can I generate the green line then? Thanks.

(ps. Any other programming languages other than MATLAB is welcome as well)

• Do you know the formula of Intensity? Implement it! – user21232 Nov 26 '16 at 6:01
• @Sardar_Usama Frankly I quite confused with the definition of intensity. In the definition en.wikipedia.org/wiki/Sound_intensity, it's defined as "the sound power per unit area", whose unit is "watt per square meter (W/m2)". But what I want is the one with unit of "dB". – walkerlala Nov 26 '16 at 6:12
• Then do the conversion! – user21232 Nov 26 '16 at 6:15
• @Sardar_Usama can you explain the conversion procedure in a formal answer? Thanks – walkerlala Nov 26 '16 at 6:21
• The wikipedia link that you mentioned explains the conversion procedure under the heading "Sound intensity level" – user21232 Nov 26 '16 at 6:30

Any audio recording is generally assumed to have units that are proportional to Pascals (Pa). This is the instantaneous pressure of the sound wave above (positive) or below (negative) the ambient pressure.

(To actually get Pascals from the WAV file, you need to know the calibration factor for your recording system.)

Loudness of sound (in air) is generally reported as Sound Pressure Level (SPL) in decibels. SPL (in dB) and pascals are related as

 Pa_ref = 20e-6; %reference pressure for SPL SPL_dB = 10*log10( (Pa / Pa_ref).^2 ); 

If you use the instantaneous pressure values (ie, the values straight from the WAV file), this equation will give you the instantaneous SPL. Based on the graph that you showed, you don't want the instantaneous SPL. You want a time-smoothed SPL.

To get a time-smoothed SPL, the standard approach is to extract the envelope of the pressure signal. There are lots of ways of doing this. One way is to rectify and low-pass filter the signal. In this case, your processing might look like this:

 %load and calibrate the data [wav,fs_Hz]=audioread(myfilename); %load the WAV file my_cal_factor = 1.0; %the value for your system to convert the WAV into Pascals wav_Pa = wav * my_cal_factor;   %extract the envelope smooth_sec = 0.125; %"FAST" SPL is 1/8th of second. "SLOW" is 1 second; smooth_Hz = 1/smooth_sec; [b,a]=butter(1,smooth_Hz/(fs_Hz/2),'low'); %design a Low-pass filter wav_env_Pa = sqrt(filter(b,a,wav_Pa.^2)); %rectify, by squaring, and low-pass filter   %compute SPL Pa_ref = 20e-6; %reference pressure for SPL in Air SPL_dB = 10.0*log10( (wav_env_Pa ./ Pa_ref).^2 ); % 10*log10 because signal is squared   %plot results figure; subplot(2,1,1); t_sec = ([1:size(wav_Pa)]-1)/fs_Hz; plot(t_sec,wav_Pa); xlabel('Time (sec)'); ylabel('Pressure (Pa)');   subplot(2,1,2) plot(t_sec,SPL_dB); xlabel('Time (sec)'); ylabel('SPL (dB)'); yl=ylim;ylim(yl(2)+[-80 0]); 

• the code is fine, but i might suggest filtering with some weighting-filter (like A-weighting) if perceived loudness, not just SPL is what is desired. coefficients for two cascaded biquads might be found here. – robert bristow-johnson Nov 27 '16 at 5:32
• I totally agree with filtering to A-weight the results. My answer was already getting pretty long, so I left it out. I totally agree, though. – chipaudette Nov 27 '16 at 23:21
• note also that you are square rooting wave_env_Pa and squaring it right away. – robert bristow-johnson Nov 28 '16 at 7:18
• I did that for clarity...especially when comparing against my initial definition of SPL_dB, which has Pascals be squared. On a PC (Matlab was the requested format), the wasted cycles of sqrt followed by square are probably inconsequential. On an embedded platform, it would indeed be stupid to do the operation as I've shown it. That's a good observation. – chipaudette Nov 28 '16 at 19:50
• thanks for your answer. Can I ask one more: What is the reason of squaring the wav_Pa? If wav_pa1 = 10 Pa and wav_pa2 = -10 Pa, then squaring both would give you 100 and 100, which are the same, but initially they have a distance of 20. After squaring them, they are now the same, but initially they are different. What is the reason of ignoring this different ? Thanks. – walkerlala Nov 29 '16 at 3:19