poctave function is providing an unexpected center frequency value and resulting amplitude
7 views (last 30 days)
Show older comments
Hi,
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
[p_oct,cf_oct]=poctave(Data(OctStartIndex:OctEndIndex),fs,'FrequencyLimits',[fBandMin,floor(fs/2)],'BandsPerOctave',bpo);
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
set(gca,'XTick',xPlot)
set(gca,'xticklabel',round(10.^get(gca,'Xtick')))
ylim([min(yPlot)-10,max(yPlot)+10])
xtickangle(90);
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')
3 Comments
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,
Colin
Answers (2)
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.
Regards,
Pratyush.
3 Comments
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)
2 Comments
ngoc quy hoang ngoc quy
on 26 Feb 2021
exactly, you can write code overrall of 1/3 octave as show below (bar color red)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/532144/image.png)
See Also
Categories
Find more on Measurements and Spatial Audio in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!