MATLAB Answers

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

68 views (last 30 days)
Matthew Casiano
Matthew Casiano on 1 Feb 2021
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
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

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.
Regards,
Pratyush.
  2 Comments
Matthew Casiano
Matthew Casiano on 24 Feb 2021
FYI, I was thinking more about how to address this. Part of the confusion might be that it is not clear whether 'FrequencyLimits' includes the band of the specified values or exlcudes the band of the specified value. I think the intent is to be inclusive since that is how it seems to behave. The caveat (which is the fix needed) would be that if the upper edge band frequency of the highest band specified exceeds Nyquist frequency, then the band should not be calculated. I think it would be helpful to note the inclusivity and caveat in the the 'FrequencyLimits' description. A more general update could include a name-value flag to specify 'inclusive' or 'exclusive' of the frequnecy limit value.

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)
  2 Comments

Sign in to comment.

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!