0
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Sorry, I thought it would be self explanatory for the fft but I haven't thought about the others. 1) Yes, the "+" operator is overloaded so the signal is multiplied with the window. 2) It's defined "20 * log10(value)" in the library. 3) I'll think about it when I get back to work, but I don't think that's a problem in my case. $\endgroup$ – Khaz42 Jun 25 '15 at 21:23
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ About the imaginary part set to 0, I've read it somewhere that it should be done that way, but the results aren't much different in both cases. Anyway, I think I should take a look at those "finite slope filters". Thanks for the feedback. $\endgroup$ – Khaz42 Jun 26 '15 at 9:57
  • $\begingroup$ The only case when just the real part of the spectrum will be sufficient to estimate the magnitude, is if your windowed time signal is even-symmetric, and can thus be decomposed into only cosines. This won't be the case when processing audio. $\endgroup$ – kippertoffee Jun 26 '15 at 10:20
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.