How to smooth FFT data on logarithmic frequency bins (octave bands)?

[Matthias Brennwald]

Hi all

I calculated the frequency response of a loudspeaker from the impulse response. 
Now I'd like to smooth the response, e.g. in a 1/24 octave band (octave, not 
Octave). I do this by defining the frequency bins (e.g., 1/24 octave bins) and 
then loop through these bins, find the data in the frequency range of each bin, 
and calculate the means of each bin (see code example below). This (i) slow and 
(ii) not elegant.

Does someone have an idea on how to improve this? Is there a way to transform 
the data such that the loop can be avoided?

In addition, it would be nice to use a more elaborate method than just 
calculating the means of each bin. For instance, something similar to a 
Savitzky-Golay smoothing would be nice.

Thanks for your help!

Matthias

=============================
Code example of my current approach:
=============================

% raw frequency response (before smoothing):
f = 10:10:40000;                        % frequency values of loudspeaker 
frequency response
Y = 10 + rand(size(f));         % frequency response (absolute values of 
Fourier transporm)

% 1/24 octave smoothing:
sm = 1/24;                                                              % 
smoothing band width (fraction of an octave)
n  = floor ( log2(f(end)/f(1)) / sm );  % number of frequency bins for avaraging
YY = repmat (NaN,1,n); ff = repmat (NaN,size(YY)); % initialize vectors that 
will contain the smoothed data
f1 = f(1);      % lower end of first frequency bin
for i = 1:n     % loop through frequency bins and calculate means of each bin
        f2 = f1*2^sm;   % upper end of frequency bin
        ff(i) = sqrt(f1*f2); % center frequency of current bin
        j = find (f >= f1 & f < f2); % find data in frequency bin f1...f2
        if any (j) % compute mean
                YY(i) = mean (Y(j));
        end
        f1 = f2; % prepare for next iteration
end

% remove missing / NaN entries:
i = find (~isnan(YY));
ff = ff(i); YY = YY(i);

plot (f,Y,'k-','linewidth',1 , ff,YY,'r-','linewidth',5);
legend ('Raw','Smoothed 1/24 octave’)

=============================