MATLAB Answers

poctave() returns incorrect fractional-octave-band results

75 views (last 30 days)
Witold Waldman
Witold Waldman on 14 Dec 2020
Answered: Pratyush Roy on 4 Feb 2021
I have encountered a bug when using the fractional-octave-band output of the poctave() function with 'psd' data. Under some circumstances, some of the bands contain significantly incorrect power levels, and sometimes two bands appear when there should only be one band present.
I have created a small M-file to demonstrate the problem, bug_poctave_psd.m, and it is attached. The output from this script is copied below. The error is easily understood when the sum of the powers in the fractional-octave-bands does not equal the total power in the input signal! the test program can easily be modified to produce errors where the results returned by poctave() are 100% out, which is clearly significant.
The error seems to occur because poctave() does not correctly handle partial FFT bin widths when performing the integration of the PSD data that falls within a fractional-octave-band. It does not correctly apportion the energy contained in a single FFT bin based on the fraction of the FFT bin that needs to be considered when a bandedge lies somewhere within the bandwidth of the requisite FFT bin. This error will be occurring for all the fractional-octave-bands, but it only becomes visible with certain test signals.
For example, here is an FFT bin, whose limits are denoted by |, with a bandedge (between two contiguous fractional-octave-bands) occurring within its confines, denoted by +.
|=====+==========|
In the line diagram above, going from left to right, it is seen that the PSD bin's contribution of energy to the left-hand fractional-octave band is less than its contribution to the right-hand fractional-octave-band. This needs to be taken into account at all upper and lower bandedges for each fractional-octave-band.
However, poctave() appears to take all the energy in the PSD bin and assign it to BOTH the upper and lower fractional-octave-bands. This is clearly incorrect, but, owing to typical FFT resolutions, the error is difficult to identify in most types of output. My simple test signal, which is a valid time-domain signal, produces significant errors in the output of poctave().
+===========================================================+
| TEST CASE 1: No cosine signal at Nyquist frequency (Fs/2) |
+===========================================================+
Number of waves = 8
Frequencies X = 125.000 250.000 500.000 1000.000 2000.000 4000.000 8000.000 16000.000
Amplitudes X = 1.414 1.414 1.414 1.414 1.414 1.414 1.414 1.414
Total Power X = 8.000
Total Power Pxx PSD = 8.000
Delta F in FFT = 1.000 Hz
Octave Fraction b = 1.000
Total Power poctave = 8.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Octave Fraction b = 1.500
Total Power poctave = 9.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
Octave Fraction b = 2.000
Total Power poctave = 10.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
Octave Fraction b = 3.000
Total Power poctave = 8.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Octave Fraction b = 6.000
Total Power poctave = 10.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
Octave Fraction b = 12.000
Total Power poctave = 10.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
Octave Fraction b = 24.000
Total Power poctave = 10.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
Octave Fraction b = 48.000
Total Power poctave = 11.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
Octave Fraction b = 96.000
Total Power poctave = 14.000
Total Power Pxx PSD = 8.000
Total Power X Data = 8.000
Incorrect power computed using results from poctave().
+=================================================================+
| TEST CASE 2: Includes cosine signal at Nyquist frequency (Fs/2) |
+=================================================================+
Number of waves = 9
Frequencies X = 125.000 250.000 500.000 1000.000 2000.000 4000.000 8000.000 16000.000 32768.000
Amplitudes X = 1.414 1.414 1.414 1.414 1.414 1.414 1.414 1.414 1.414
Total Power X = 10.000
Total Power Pxx PSD = 10.000
Delta F in FFT = 1.000 Hz
Octave Fraction b = 1.000
Total Power poctave = 8.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Incorrect power computed using results from poctave().
Octave Fraction b = 1.500
Total Power poctave = 9.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Incorrect power computed using results from poctave().
Octave Fraction b = 2.000
Total Power poctave = 10.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Octave Fraction b = 3.000
Total Power poctave = 8.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Incorrect power computed using results from poctave().
Octave Fraction b = 6.000
Total Power poctave = 10.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Octave Fraction b = 12.000
Total Power poctave = 10.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Octave Fraction b = 24.000
Total Power poctave = 10.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Octave Fraction b = 48.000
Total Power poctave = 11.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Incorrect power computed using results from poctave().
Octave Fraction b = 96.000
Total Power poctave = 14.000
Total Power Pxx PSD = 10.000
Total Power X Data = 10.000
Incorrect power computed using results from poctave().
>>
  5 Comments
Mathieu NOE
Mathieu NOE on 17 Dec 2020
hello
maybe you should send your test case to the TMW for double check
all the best

Sign in to comment.

Answers (1)

Pratyush Roy
Pratyush Roy on 4 Feb 2021
Hi Witold,
I have heard that this issue is known and the concerned parties may be investigating further.
Regards,
Pratyush.

Community Treasure Hunt

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

Start Hunting!