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

1 次查看（过去 30 天）
Matthew Casiano
Matthew Casiano 2021-2-1
评论： Matthew Casiano 2021-12-3
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 个评论
显示 1更早的评论
Matthew Casiano
Matthew Casiano 2021-2-1
Thanks for your comment. I agree it seems like a bug. However, I don't necessarily want to specify a guard band such as 2.56x. In general, I'd like to keep the capability of a band reaching the Nyquist frequency if it can, but not using a partial band. In addition to processing simulation data, our antialiasing filters have a very good response near the Nyquist frequency. So in general it would be desired if the octave band can reach the Nyquist frequency if it is within range. I'm guessing to account for the bug though, an appropriate max frequency limit would have to be defined - similar to your approach, but maybe based on the octave band definition logic.
Colin Brown
Colin Brown 2021-2-8
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Pratyush Roy
Pratyush Roy 2021-2-8
Hi Matthew,
I have brought this issue to the notice of our developers. They will investigate the matter further.
Regards,
Pratyush.
  3 个评论
显示 1更早的评论
Matthew Casiano
Matthew Casiano 2021-2-24
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.
Matthew Casiano
Matthew Casiano 2021-12-3
@Pratyush Roy
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,

请先登录，再进行评论。

ngoc quy hoang ngoc quy
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 个评论
Matthew Casiano
Matthew Casiano 2021-2-25
It will depend on the units of your data. My data are in the EE Units system, so I have pref first converted to psi. Use 2e-5 Pa if you're in SI. But for convention, my plot shows units as dB re 2e-5 Pa.
I guess for clarity the commented line I have in the above script for the data vector should show psi rather than engineering units: % data vector (psi)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Octave 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by