# ConstantQ transformation method

I've got this transform method with a triple-nested loop. The Generate methods do their own caching so they are fast. K will have a worst-case value of 150. N will have a worst-case value of about 20000. It's that 200002 that's killing me. It looks to me like the innermost loop could be moved out, but would you store every possible sum?

    public static void Transform(IList<float> x, int offset, int sampleRate, int binsPerOctave,
double minFrequency, double maxFrequency, out double[] magnitudes, out double maxMagnitude)
{
var Q = 1.0 / (Math.Pow(2.0, 1.0 / binsPerOctave) - 1.0);
var K = (int)Math.Ceiling(Math.Log(maxFrequency / minFrequency, 2.0) * binsPerOctave);

magnitudes = new double[K];
var maxN = (int)Math.Round(Q * sampleRate / minFrequency);
var halfN = -maxN / 2;
double[] wcs = GenerateWindowConstants(maxN);
maxMagnitude = double.NegativeInfinity;

for(int k = 0; k < K; k++)
{
var N = (int)Math.Round(Q * sampleRate / (minFrequency * Math.Pow(2.0, (double)k / binsPerOctave)));
Complex[] piqs = GeneratePiQConstants(N, Q);

var maxMag = double.NegativeInfinity;
for(int i = 0; i < N; i++)
{
Complex sum = new Complex();
for(int n = 0; n < maxN; n++)
{
var index = n + offset + halfN;
index = Math.Min(x.Count - 1, index);
index = Math.Max(0, index);
sum += (wcs[n] * x[index]) * piqs[(n + i) % N];
}
sum /= maxN;
maxMag = Math.Max(maxMag, sum.Magnitude);
}
magnitudes[k] = maxMag;
maxMagnitude = Math.Max(maxMagnitude, maxMag);
}
}

• What is the aim of this code?
– Oli Charlesworth
Apr 5, 2012 at 19:41
• This code is a modified ConstantQ transform -- modified in a way to try all possible phase shifts and use all the data available for higher frequencies.
– Brannon
Apr 5, 2012 at 19:49
• How much speedup do you look for? What is the hot path? Have you tried libraries written in C++ that take advantage of platform? Apr 5, 2012 at 19:57
• I'm looking to speed this thing up from hours to minutes. The hotpath is the sum += ... I've thought about trying out the C++ compiler to see if it would give me any significant help. I've also thought about ditching the Complex struct and going with my own math as I know the CLR does a lot better optimization with simple types.
– Brannon
Apr 5, 2012 at 20:31

It looks like based on the 2 inner loops where you are calculating maxMag you could refactor the code to utilize Parallel.ForEach() so you can start to use multiple cores for the math. Just put all of the individual maxMags into a concurrent dictionary. When your done, run your aggregation methods. You'll have to decide whether it is better as an array that gets copied, or whether to refactor them into a concurrent collection type. I hope this helps.

Edit/Update:

Per your comment relating to using all cores "However, the best improvement this would bring on typical hardware would be 8x. 8x isn't quite going to be enough for my application." (found under cgilmeanu's answer)

If utilizing all of the computer resources and taking an 8-fold gain won't get anywhere close, then ...

Answer = you need to horizontally scale across multiple computers until you reach the point where you are satisfied with the "timed" outcome.

Going to c, c++, or .net-under-the-hood should help make gains for you. However, commonly, the programming involved results in something complicated, hard to understand, and more extensive documentation. This can make it tough to be nimble and make changes. Also, if you can assume that the computational load/iterations may increase in the future, your extensive efforts in extremely fine tuning this item may be all for nothing since it was internally just vertical scaling by optimization.

Note the comment on FFT here.

(I assume your hot path is sum += (wcs[n] * x[index]) * piqs[(n + i) % N];).

If N is a power of 2, you could use a shift instead of a divide. Have you tried padding up to a power of 2?

Also, are you following the advice in this paper?

• I definitely need to study that paper more. I'm aware that ConstantQ can be done with FFT. However, ConstantQ differs from my code in that it only has the one nested loop. (eg. replace the 'n' in my code with 'i'.) I haven't dug through the math to see how my current code could utilize the FFT (and I'm not sure my math skills are up to that task).
– Brannon
Apr 5, 2012 at 20:29

I rewrote the inner most loop using parallel functions. I don't have all necessary classes for running tests on my own, so if you want to try it just let me know the results. Hope this helps.

You could replace the following code

Complex sum = new Complex();
for(int n = 0; n < maxN; n++)
{
var index = n + offset + halfN;
index = Math.Min(x.Count - 1, index);
index = Math.Max(0, index);
sum += (wcs[n] * x[index]) * piqs[(n + i) % N];
}


with this one

//ReaderWriterLockSlim rwlock = new ReaderWriterLockSlim();
Complex sum = new Complex();
object lockObj = new object(); // this can be put outside the loop, in order not to create it every time the index advance
ParallelOptions options = new ParallelOptions { MaxDegreeOfParallelism = 4 }; // this can be put outside the loop, in order not to create it every time the index advance
Parallel.For(0,
maxN,
options,
() => 0,
(int n, ParallelLoopState loopState, Complex threadLocalSum) => {
var index = n + offset + halfN;
index = Math.Min(x.Count - 1, index);
index = Math.Max(0, index);
threadLocalSum += (wcs[n] * x[index]) * piqs[(n + i) % N]);