Unexpected results using octaveFilter in Matlab

Question

I'm aiming at obtaining a 1/3 octave band frequency signal for specific centre frequencies. I want to use Matlab's octaveFilter functionality but I was expecting to get a single 1/3 octave band peak at 1000 Hz, but instead far left and far right of 1000 Hz very positive sound pressure levels were computed.

What do I do wrong?

    Fs = 48000;                         % Sampling rate
    t = 0:1/Fs:1-1/Fs;                  % Time vector of 1 second
    f = 1000;                           % Frequency of signal
    dpres = sin(2*pi*f*t);              % Signal in [Pa]

    frCen = [100 300 600 800 1000 1200 1300 1600];

    [Spl,frCen] = CompOctSplFreq(dpres,frCen)

    figure()
    semilogx(frCen,Spl,'ro-')
    title('1/3-Octave Filtered SPL over Frequency')
    xlabel('Center Frequency of Octave Band Filter [Hz]')
    ylabel('SPL [dB]')        

            function [Spl,frCen] = CompOctSplFreq(dpres,frCen);

                % Reads singal in pressure over time
                freqNum = length(frCen);
                Fs = length(dpres);

                % Setting up the 1/3-octave filters for all center frequencies
                for f = 1:freqNum
                    octaveFilterBank{f} = octaveFilter(frCen(f),'1/3 octave','SampleRate',Fs,'FilterOrder',12);
                    % Filtering the signal with the corresponding filters
                    dpresFiltered(:,f) = octaveFilterBank{f}(dpres');
                    % Getting the average for each filter frequency band
                    drms(f) = sqrt(sum(dpresFiltered(:,f).^2)/length(dpresFiltered(:,f)));
                end

                % Converting the root mean square pressure to SPL            
                pRef = 20e-06;                      % Reference pressure            
                Spl = 20*log10(drms/pRef); 


            end

Answer 1

The main issue is that you are computing the power of signals which in some cases include a significant transient response. To illustrate this the following will go through the derivation of the expected steady-state response, and compare the resulting computed responses including and excluding the transient portion.

Based on the algorithm description of octaveFilter, it is possible to construct the ideal steady state response of each filter in your filter bank. These responses are given in the following graph:

Note that the responses do go down significantly outside their specified band, but they do not completely vanish. As a result input signals with frequency outside the specified band will still cause some output (albeit small) on all the filters in the filter bank. For a sinusoidal input signal at 1000Hz, you could get the steady state response at the output of each filter by looking up the intercept of the corresponding filter response with the 1000Hz vertical line in the above graph. This should give you the following graph for the response of the filter bank to a 95dB input at 1000Hz:

Now you may notice that your results are pretty close to these steady state responses for the filters centered near 1000Hz, but the difference increases as we get further away from 1000Hz. If we look at one of the filtered signal in the time-domain, we may notice that indeed a significant portion of energy can be found during the initial transient of the signal:

Simply getting rid of the transient should bring you back closer to the expected steady-state response. The trickier part is to figure out how many samples need to be discarded. As a general rule of thumb I generally use approximately 5 time constants, where one time constant is on the order of the inverse of the passband filter half bandwidth. This can be computed with the following:

G = 10^(3/10);
b = 3; % for 1/3-octave
fhi  = frCen(f)*G^(0.5/b);
flow = frCen(f)*G^(-0.5/b);
M = 10/(fhi-flow); % 5 time-constants (one time-constant is ~ 1/((fhi-flow)/2) )
drms(f) = sqrt(sum(dpresFiltered(M:end,f).^2)/length(dpresFiltered(M:end,f)));

Comparing the resulting response excluding the initial samples corresponding to the transient with your previous result results would show that this indeed had a significant impact on the graph, and that the results are now in much better agreement with our earlier ideal steady-state response graph: