I see you're still falling into the abyss of C++ metaprogramming. I still insist that using metaprogramming in signal processing is a bad idea. But I think I'll give some advices to this example. First of all, move out from the cycle set of constant calculations:

Code: Select all

`template <typename CalcType>`

CalcType BandlimitedSawtoothWave<CalcType>::operator()()

{

CalcType sample = {};

CalcType k = CalcType{2.0} * M_PI * phase.GetValue();

for (std::size_t i = 1; i <= numharmonics; ++i)

{

sample += std::pow(CalcType{-1.0}, i) * std::sin(k * i) / i;

}

sample = CalcType{-1.0} / M_PI * sample;

phase += phasedelta;

return sample;

}

Now we see pretty expensive operation: std::pow that uses -1 as an argument.

Let's evaluate it's values:

std::pow(-1, 1) = -1

std::pow(-1, 2) = 1

std::pow(-1, 3) = -1

std::pow(-1, 4) = 1

Oh damn, we can use just multiplier -1 for odd i and +1 for even i:

Code: Select all

`template <typename CalcType>`

CalcType BandlimitedSawtoothWave<CalcType>::operator()()

{

CalcType sample = {};

CalcType k = CalcType{2.0} * M_PI * phase.GetValue();

for (std::size_t i = 1; i <= numharmonics; ++i)

{

if (i & 1)

sample -= std::sin(k * i) / i;

else

sample += std::sin(k * i) / i;

}

sample = CalcType{-1.0} / M_PI * sample;

phase += phasedelta;

return sample;

}

Now we see that we can eliminate the (i & 1) condition by unrolling the cycle:

Code: Select all

`template <typename CalcType>`

CalcType BandlimitedSawtoothWave<CalcType>::operator()()

{

CalcType sample = {};

CalcType k = CalcType{2.0} * M_PI * phase.GetValue();

std::size_t i=1;

size_t loops = numharmonics - 1;

while (i < loops)

{

sample -= std::sin(k * i) / i;

i++;

sample += std::sin(k * i) / i;

i++;

}

if (numharmonics & 1) // numharmonics is odd

sample -= std::sin(k * i) / i;

sample = CalcType{-1.0} / M_PI * sample;

phase += phasedelta;

return sample;

}

Then we can calculate sine and cosine once before cycle and then use complex vector rotation. That works faster because requires less math.

We can calculate the initial complex vector for i=1:

Code: Select all

`double v_re = std::cos(k);`

double v_im = std::sin(k);

And the rotation vector will be the same because the rotation angle is k:

Code: Select all

`double r_re = v_re;`

double r_im = v_im;

Now our code becomes looking like this:

Code: Select all

`template <typename CalcType>`

CalcType BandlimitedSawtoothWave<CalcType>::operator()()

{

CalcType sample = {};

CalcType k = CalcType{2.0} * M_PI * phase.GetValue();

double v_re = std::cos(k), v_im = std::sin(k);

double r_re = v_re, r_im = v_im;

double t_re, t_im; // Temporary vector

std::size_t i=1;

size_t loops = numharmonics - 1;

while (i < loops)

{

sample -= v_im / i; // Using the imaginary part of the R vector because it's sine

// Rotate the vector { v_re, v_im } by multiplying on {r_re, r_im} and store into { t_re, t_im }

t_re = v_re * r_re - v_im * r_im;

t_im = v_re * r_im + v_im * r_re;

sample += t_im / (i + 1); // Using the imaginary part of the T vector because it's sine

// Rotate the vector { t_re, t_im } by multiplying on {r_re, r_im} and store into { v_re, v_im }

v_re = t_re * r_re - t_im * r_im;

v_im = t_re * r_im + t_im * r_re;

i += 2;

}

if (numharmonics & 1) // numharmonics is odd

sample -= v_im / i; // We don't need to rotate the vector because it's the last calculation

sample = CalcType{-1.0} / M_PI * sample;

phase += phasedelta;

return sample;

}