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;
}