poctave function is providing an unexpected center frequency value and resulting amplitude
조회 수: 8 (최근 30일)
이전 댓글 표시
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
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
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
ngoc quy hoang ngoc quy
2021년 2월 25일
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
ngoc quy hoang ngoc quy
2021년 2월 26일
exactly, you can write code overrall of 1/3 octave as show below (bar color red)
참고 항목
카테고리
Help Center 및 File Exchange에서 Octave에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!