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

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.

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.

• 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. – Khaz42 Jun 25 '15 at 21:23

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.

• 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. – Khaz42 Jun 26 '15 at 9:57
• 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. – kippertoffee Jun 26 '15 at 10:20

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)