EQ - Calculate the power per Octave of a signal through FFT

Question

I'm working on a audio Equalizer, and I'm not sure if the way I obtain the power per octave of my signal is good.

Here is the code I have :
(I use the third octave bands, which is initialized earlier in the program)

void    Record::ProcessData(short* Buffer, size_t BufferSize)
{
    /* BufferSize is usualy equal to 131072 */
    /* Initialize the signal and the Hanning Window */
    Aquila::SignalSource signal = Aquila::SignalSource(Buffer, BufferSize, 44100);
    Aquila::HannWindow HannWindow(BufferSize);

    /* FFT */
    /* The library overload the + operator, so the signal is multiplied with the window */
    auto fft = Aquila::FftFactory::getFft(BufferSize);
    Aquila::SpectrumType spectrum = fft->fft((signal + HannWindow).toArray());

    /* Set the imaginary part to 0 */
    for (size_t i = 0; i < BufferSize; i++)
        spectrum[i]._Val[1] = 0;

    /* OctaveRange is a class with the minimum, center and maximum frequency of an octave band */
    list<OctaveRange>::iterator it;
    map<double, double> sumMap;
    double db, f, key;

    /* Initialize the map values to 0 */
    for (it = thirdOctave.begin(); it != thirdOctave.end(); it++)
        sumMap[(*it).getCtr()] = 0.0;

    /* The calcul of the power per octave starts from here */
    for (size_t i = 0; i < (BufferSize / 2); i++)
    {
        /* f = frequency at bin i */
        f = (double(i * 44100) / BufferSize);
        for (it = thirdOctave.begin(); it != thirdOctave.end(); it++)
        {
            key = (*it).getCtr();
            if (f >= (*it).getMin() && f < (*it).getMax())
                sumMap[key] += pow(spectrum[i].real(), 2) + pow(spectrum[i].imag(), 2);
        }
    }

    for (it = thirdOctave.begin(); it != thirdOctave.end(); it++)
    {
        key = (*it).getCtr();
        /* db is the power (in dB) per third octave. */
        db = 20.0 * std::log10(sqrt(sumMap[key]);
    }
}

When I record pink noise with this I obtain some believable values, but I'm still not sure if this is because the code is good, or if it's a happy coincidence.

So my question is simple :
Is this the good way to obtain the power per octave band of a signal ?

Feel free to ask if something is not clear enough.
Thanks for your time.

Answer 1

Difficult to say without a reference to the library functions and classes that you are using. In general it looks okay (if fairly inefficient).

A few things to look at:

1. The Hamming window needs to be multiplied with the signal, not added to it. Whether the code is correct really depends on how the "+" operator for the clasess is defined
2. db() can be 10*log or 20*log, whether you need the sqrt depends again on how its's defined in the library
3. Your octave band filters are bascially "brick wall" filters. Most spectrum analyzers do have a finite slope filters so any one frequency contributes to multiple bands (just by different amounts). Whether that's a problem or not depends on the application

Re point 3: consider the output of the analyzer when you input a sine wave and slowly sweep the frequency from 1kHz to 2 kHz. With a finite slope filter you will see the energy in the 1 kHz octave band slowly going down and the one in the 2 kHz band slowly coming up. They will be equal when the frequency is 1414 Hz. In your case the 1 kHz octave band energy will be constant and once you reach 1414 Hz all the energy will immediately flip to the 2 kHz band.

Answer 2

You are setting the imaginary part of your spectrum to 0, so when you calculate the power here:

 sumMap[key] += pow(spectrum[i].real(), 2) + pow(spectrum[i].imag(), 2);

it's just the squared real part. I don't see any reason for this in your application.

Also, you could save yourself a sqrt() in each bin and just do:

db = 10.0 * std::log10(sumMap[key]);

Those points aside, the general approach is sound. Parseval's theorum backs this up. Although, as Hilmar stated, usually fractional octave band filters in analysers are not brick wall filters.

Answer 3

If setting the imaginary part to zero does not make much difference, it's coincidental.

Here's a page with some discussion of how to sum the FFT bins for octave analysis:

http://www.mstarlabs.com/docs/tn257.html

A more sophisticated way to sum the FFT bins to get the octave band amplitudes is to use the constant Q transform. You can find references on the Wikipedia page:

https://en.wikipedia.org/wiki/Constant_Q_transform

Reference 3 has the math.

In the "traditional" method you apply a filter bank for the top octave, decimate by 2, apply the same filter bank to the top octave of the result, and so on. For a MATLAB implementation and references see (PDF)

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.57.5728&rep=rep1&type=pdf