poctave function is providing an unexpected center frequency value and resulting amplitude

3 views (last 30 days)
I am looking for feedback on the poctave function. I have found that in some cases, the 'FrequencyLimits' flag does not hold true. I am applying the poctave function to 1/3 octaves ('BandsPerOctave' is 3). For example, for a sampling rate of 10,000 sps, I'd expect the highest 1/3 octave band would have a center frequeny of ~3981 Hz (since the upper end of that band only reachs ~4467 Hz) and since the next band centered at ~5012 Hz exceeds the Nyquist frequency. However, the output provided by the function lists values for the 5012 Hz band. Further, plotting the amplitude of this center frequency does not follow the expected trend. It appears that only a partial band is being evaluated. I wouldn't expect that a partial band would be evaluated given the FrequencyLimit since the overall value for that band would change if the rest of the band was included.
The excerpt of my example here evaluates the function and plots 1/3 octave bands for white noise. There is some variation due to the random signal input, but the last band at ~5012 Hz is present and is always lower then the expected trend.
%% Inputs
t1=1; % start time for octave band spectrum (sec)
t2=19; % end time for octave band spectrum (sec)
fs=10000; % sampling rate (sps)
fBandMin=20; % minimum low-frequency band (Hz)
bpo=3; % bands per octave, e.g., 3 is one-third octave bands
pref=2e-5/(4.4482216152605/0.0254^2); % decibel reference value (psi), 1 psi = 4.4482216152605/0.0254^2 Pa [exact]
%% Example data
TimeData=0:1/fs:20; % time vector (s)
Data=rand(1,length(TimeData)); % data vector (EU)
tStart=TimeData(1); % data start time
%% Octave Spectrum
OctStartIndex=round((t1-tStart)*fs)+1; % find index closest to specified start time
OctEndIndex=round((t2-tStart)*fs)+1; % find index closest to specified end time
% Note this is the average power over the band
spl=10*log10(p_oct/pref^2); % sound pressure level (dB re pref)
%% Plotting
% First evaluate the log10 of the center frequency values and then rewrite
% tickmarks using the base 10 antilog. This forces Matlab to produces
% equal-sized bars and then corrects the labels
xPlot=log10(cf_oct); % base 10 logarithm of the center frequencies
yPlot=spl; % spl values
bar(xPlot,yPlot) % octave band spectrum plot
grid on
xlabel(gca, 'Band Center Frequency (Hz)')
ylabel(gca, 'SPL (dB re 20\muPa)')
title(gca,{'1/3 Octave Band Spectrum'},'fontsize',14,'fontweight','bold')
Colin Brown
Colin Brown on 8 Feb 2021
I wonder why poctave also imposes a minimum frequency limit of 3 Hz. I have acoustic signals from sources with frequencies as low as 0.05 Hz which I can see using power spectra density estimates. Can this be fixed as well?
Thank you,

Sign in to comment.

Answers (2)

Pratyush Roy
Pratyush Roy on 8 Feb 2021
Hi Matthew,
I have brought this issue to the notice of our developers. They will investigate the matter further.
Matthew Casiano
Matthew Casiano on 3 Dec 2021
Hi, I received a noticed from Mathworks Technical Support that this reported issue has been addressed in the R2021b release. I see some other similar poctave issues have been resolved in this release, but this one does not appear to have been resolved. Running the script above still shows the same issue as before (evident in the last octave band after running the script).
In the example above, entering the upper 'FrequencyLimits' value as the Nyquist frequency (sampling rate = 10,000 sps; Nyquist = 5,000 Hz) causes poctave to produce a band at 5012 Hz. Some type of 'octave smoothing' (similar to what is described in the documentation for frequency bins) would be needed to account for the fraction of power corresponding to the percentage of the data available within the band. For this specific case, the last 1/3rd octave center frequency at 5012 Hz has a lower and upper edge at 4467 Hz and 5623 Hz; but the data is only available to Nyquist at 5000 Hz - so some type of octave smoothing is necessary if this band is included.
Thank you,

Sign in to comment.

ngoc quy hoang ngoc quy
ngoc quy hoang ngoc quy on 25 Feb 2021
follow me, spl=10*log10(p_oct/pref^2) with pref = 2*10^-5 (pa) don't pref=2e-5/(4.4482216152605/0.0254^2) (psi)

Sign in to comment.




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by