poctave() return value for acoustics analysis

Hi all,

I need some assistance understanding the "units" of the data that is returned when calling the MATLAB function [p,cf] = poctave();

I am analysing an audio signal that has been recorded with a microphone. The data is imported into the MATLAB workspace and converted to a calibrated pressure value with units of Pascals (Pa). So the units of the data that I pass to poctave() are Pascals (Pa). The following code snippet shows how I'm using poctave().

flims = [20 Fs/2];	% set the frequency limits of my analysis
bpo = 3;		% Third-Octave analysis required

opts = {'FrequencyLimits',flims,'BandsPerOctave',bpo};	% prepare these options for passing to poctave()

% apply the 3rd-Octave filter-bank to my data
[p, cf] = poctave(pressureData,Fs,opts{:});

(Where "pressureData" is my audio data to be filtered, and "Fs" is the sampling frequency of the data.)

My first question is:

1) What are the units of the octave spectrum data returned unto variable 'p'?

I need to know this for my conversion of the data into a dB value using 10*log10(p/pref).

I understand that if I use poctave() as follows:

 poctave(pressureData,Fs,opts{:}) 

then by default it assigns "pref" a value of 1 and the results are as in the following image:

Without units though, the result seems meaningless...

My second question is:

2) what is pref in the previous equation to ensure correct dB values that correlate with a sound pressure level? do I use:

pref = 0.00002;     % 20 micro Pascals reference sound pressure level

or:

pref = 0.00002^2;     % 20 micro Pascals (squared); i.e. is p returned from poctave() a pressure squared value or a power value??

I want to create a resulting 3rd-Octave spectrum whos dB values match those as though it was being displayed on a Sound Level Meter. I will then correlate these values with those on our Sound Level Meter to validate the data to be processed within MATLAB.

Thanks so much for any assistance you can offer!!

[Mickey Yoseph]

Hello mate - did you get an answer to this in the end? In the exact same boat as you for both questions.

Require dB SPL but clearly poctave returns negative dB values which aren't helpful so need 'pref' to be 0.0002.

Any help?

[zein]

i think that the output in the graph is 10log(p) where p is the spectrum in 1/3 octave and to calculate the SPL-1/3 OCTAVE

SPL-1/3octave=10log(p/pref.^2)

[Michiel Heyns]

Because poctave calculates the power over the octave frequency limits from the PSD, it already quadrates the signal. Therefore, to see Sound Pressure Level on the graphs created by poctave, the inputs should be:

x=Pressure/P0;

[p,f]=poctave(x,Fs,...)

Note, the algorithm calculates the output as follows:

p=10^(S/10)

Therefore, the output can be interpreted as sound power or sound intensity, not pressure.

[Michiel Heyns]

To then calculate the leves:

LI=10*log10(I/I0)

[Michiel Heyns]

You can easily test the algorithm with a sine wave. See below:

% This program test the accuracy of the sounanalysismain program with known
% signals.
% Programmed by: Dr Michiel Heyns
% 2024-08-24

% Functions used
LI=@(Lp) 10*log10(sum(10.^(Lp/10)));  % To calculate the equivalent sound pressure level.

% Define parameters and frequencies.
P0=20e-6;  % The reference pressure [Pa].
I0=1e-12; % The reference sound intensity [W/m^2].
Fs=98304; % Sampling frequency of the signal [Ha].
F=500;  % Frequency of pressure signal.
dB=90; % The dB value to be constructed.  Note that rms values are used in L_eq calculations.

ti=0;te=5; t=linspace(ti,te,(te-ti)*Fs);
P=10.^(dB/20)*P0;  % But, this is rms amplitude.  To determine the amplitude of the sine wave:
P=sqrt(2)*P; % The amplitude of the sine wave that will give the rms pressure above.
Pressure=P*sin(2*pi*F*t);

X=Pressure/P0;
options1={'FrequencyLimits',[30 16000],'BandsPerOctave',3,'Weighting','A','FilterOrder',8};
options2={'FrequencyLimits',[30 16000],'BandsPerOctave',3,'FilterOrder',8};
options3={'spectrogram','FrequencyLimits',[30 16000],'BandsPerOctave',48,...
    'OverlapPercent',50,'Weighting','A','FilterOrder',8};

[Lp1,f] = poctave(X,Fs,options1{:});
[Lp2,f] = poctave(X,Fs,options2{:});
[Lp3,f] = poctave(X,Fs,options3{:});
Leq1=round(LI(10*log10(Lp1)),1);
Leq2=round(LI(10*log10(Lp2)),1);
Leq3=round(LI(10*log10(Lp3)),1);
% To see the octave plot, use the command as follows:
figure
poctave(X,Fs,options1{:})
title('1/3 - Octave Spectrum, A-weighted, L_eq = ', num2str(Leq1))
figure
poctave(X,Fs,options2{:})
title('1/3 - Octave Spectrum, No weighting, L_eq = ', num2str(Leq2))
figure
poctave(X,Fs,options3{:})
title('1/48 - Spectrogram, A-weighted, L_eq = ', num2str(Leq3))

Answer 1

[Mark Thompson]

OK. Thanks for the reminder of this one @Mickey Yoseph. I'll post the solution I ended up deriving.

So, as I mentioned in my original question I had an array of calibrated pressure values. (My final solution included the creation of a structure called "source" to store all relevant information about my source file, hence why you see it referenced in the following):

pref = 0.00002;     % reference pressure in air (Pa)

% Setup octave-band filter parameters:
flims = [20 source.fs/2];   % in this case 20 Hz to 24 kHz
bpo = 3;                    % Third-Octave Filtering
opts = {'FrequencyLimits',flims,'BandsPerOctave',bpo};
% filter my source pressure data:
[source.p3rd, source.cf3rd] = poctave(source.pressure,source.fs,opts{:});
% convert filtered data to calibrated SPL values referenced to 20 uPa
source.ThirdOctDB = 10*log10(source.p3rd/(pref^2));

So, answering Question 1:

- Use poctave and get filtered pressure data in whatever Octave-band spacing you want. Note that what is returned is a pressure squared value (i.e. the average power for that Octave-band).

And answering Question 2:

- understanding that poctave returns pressure squared values we then convert to SPL figures (dB re: 20 uPa) using the equation:

Where,

I hope this helps clarify how you can use poctave and obtain referenced SPL values for acoustic pressure signals.

[Mickey Yoseph]

Thanks very much Mark this has clarified things perfectly! You're a hero!

[zein]

is it possible to change the standard octave band limits to the following bands (50-100,100-200,200-400,400-800,800-1600,1600-3200,3200-6400) with center frequency (75,150,300,600,1200,2400,4800)?

[Slimane Ouakka]

why in the code are you do not square "source.p3rd"?

[Mark Thompson]

@Slimane Ouakka, poctve() returns the pressure^2 value already, which is why I don't square it (again) in my code.

Answer 2

[ngoc quy hoang ngoc quy]

I think you are confused average power with pressure squared

[Mark Thompson]

Isn't Power proportional to the square of the pressure amplitude?

Answer 3

[andytodd.msp]

Hi All, I'm just curious as to the solution to this - as we are in a similar position (using calibrated input data and trying to have results displayed in typical dB SPL values).

I've tryed using the methods above, but I'm not sure if I'm missing something - as we do not seem to get reasonable answers using this method.

For example, the output of poctave at 1kHz is -14.7dB (Average Power). We know the SPL at 1kHz is actually 79.4dB (SPL). I'm struggling to find the link between these two values so that we can convert between the two.

Any help is much appreciated!

Andy

[Mark Thompson]

Hi Andy, are you converting your data samples (imported into MATLAB) to the units you are working in first (e.g. convert each sample to a pressure in Pascals)?

Also, I'm assuming you have a microphone calibrator sitting on top of a reference microphone and are recording that signal, then running it through the poctave function?! I assume your -14.7dB result at 1 kHz is reflecting the spread of energy due to the the windowing in the poctave(), and that is a level below 0dB. It also seems like you're missing a step of converting to actual pressure values first.

Can you confirm these things?

[andytodd.msp]

Hi Mark,

No - the data we are importing is already in Pascal units. Its coming from a calibrated measurement system.

What I'm trying to do, is have poctave display the octave values as dB re 20 micropascals - rather than the default 'dB Average Power'.

Does that make sense?

[Mark Thompson]

Hi @andytodd.msp, that's exactly what source.ThirdOctDB is when my code returns. Again, you won't get exactly 79.4 dB re 20 micropascals from this function; refer to my comment regarding windowing and the spread of energy over the width of the octave band you are looking at (hence the term "average" I guess).

Answer 4

[Aitor]

Thank you @Mark Thompson. Now the calibrated values are stored in source.ThirdOctdB. But how do you plot them?Bar plot?

Thanks in advance